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Lt is perfechiy acrue) asmthespnilosophers say, that 
life must be understood backwards. But they forget the 
other proposition, that it must be lived forwards. And 
if one thinks over that proposition it becomes more and 
more evident that life can never be understood in time 
Simply because at no particular moment can one find the 
necessary resting place from which to understand it 


—-backwards. 


-Kirkegaard, Journals, 1843 
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ABSTRACT 


This work is concerned with the general problem of 
a charged particle in a potential field of two Coulomb 
Centers. 

It is common knowledge that the Schroedinger equation 
for such a system within the Born-Oppenheimer approximation 
can be separated into a system of three simultaneous second 
order linear differential equations, if use is made of 
spheroidal coordinates. This separation of variables gives 
rise to three separation constants. Since the wave function 
has to be single valued one of these three constants is 
required to be integer. For an acceptable solution of the 
system the remaining two separation constants have to be 
equal. Furthermore they have to be eigenvalues of a pair of 
linear differential equations such that their solutions, 
commonly referred to as "generalized spheroidal wave 
functions", exist and remain bounded in the physically 
meaningful intervals of the corresponding variables. 

This thesis is divided into five chapters which 
concentrate on the following topics : (1) Review of 
existing algorithms for the numerical solution of the 


system of simultaneous linear differential equations, 
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originating mainly from Hylleraas and Jaffé. It was found 
that the radial expansion proposed by Jaffé provides no 
improvement over the previous Hylleraas expansion in 
respect to convergence properties if an equal number of 
terms is considered in both cases. 

(II) New expansions are proposed. An angular expansion 
in ultraspherical polynomials seems sufficiently rapidly 
convergent for most practical purposes. In addition an 
angular expansion in confluent hypergeometric functions 
was found which could possibly be useful if a calculation 
for a heteronuclear species with large internuclear 
distance is to be done. For the radial expansion no 
improvement of the convergence rate could be achieved. 

(III) Integral transform techniques are explored. The 
Laplace transform provides a reasonable explanation for 
the equivalence of the convergence rates of the two 
existing radial expansions, but no cure for convergence 
problems with the radial expansions could be found with 
the help of integral transforms. 

(IV) The system consisting of a particle in the dipole 
field of two oppositely charged point charges serves as 
testing ground for the algorithms available up to that 
point. The need for further improvement of the radial 
expansions is made apparent. The binding energy of a 


particle in a dipole field becomes zero at.some .critical 
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internuclear distances depending on the state considered. 


A theorem is formulated stating 
distances depend on the angular 
and it furthermore allows us to 
distances by a simple algorithm 


a root search. Numerical values 


that these critical 
quantizations exclusively 
calculate these critical 
consisting essentially of 


for ten critical 


internuclear distances are provided in the text. The 


critical distance for the ground state published 


previously by Turner and Fox is 


in agreement with the 


one calculated by using the present approach. 


(V) A radial expansion with the convergence rate 


accelerating as the "united atom" limit is approached 


is found. 
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CHAPTER I 


Introduction and Summary of Previous Work 


The theoretical investigation of the hydrogen 
molecule ion and related species is a problem which 
arises in many contemporary quantum mechanical 
treatments as a means of testing theories and 
algorithms at the ab initio level. The Schroedinger 
equation for the one-particle two-center problem in 
the Born-Oppenheimer approximation can be separated 
into a system of three simultaneous linear DE's 
(differential equations) if use is made of spheroidal 
coordinates. Two of these three DE's are closely 
related to the DE of the spheroidal wave functions 
frequently encountered in potential theory. In the 
general case these two DE's are called GSW-DE's 
("generalized spheroidal wave DE's") a notation 
widely accepted in literature. 

The algorithms for the solution of the system of 
DE's were developed at the dawn of modern quantum 
mechanics mainly by Hylleraas and galt ewmthearadial 
expansions proposed by Hylleraas and Jaffé are in this 


text referred to as "Hylleraas-" and "Jaffé-expansion". 
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They will be shown to generate algorithms with 
identical rates of convergence. This was not known 
previously. 

The actual DE of the spheroidal wave functions 
occurs if a homonuclear case or a "particle ina 
dipole field" is treated. A considerable amount of 
iInteresteis paideto a system consisting of¥a particle 
in the field of two charges with equal magnitude and 
opposite Sign i.e. a dipole field of two point 
charges. The interest in this special case is partially 
justified by the fact that from the potential curve for 
one dipole the potential curve for a different dipole 
corresponding to different magnitude of the point 
charges is obtained by simply scaling the internuclear 
distance. 

Since the separation of variables for the one-particle 
two-center problem in spheroidal coordinates is a common 
COpic in Introductory uantumechemistry courses, no 
effort has been made here to show in detail how this 
system of linear DE's is obtained from the Schroedinger 
equation in spheroidal coordinates by separation of 
variables. Instead existing algorithms for the numerical 
solution of this system of DE's are presented in this 
chapter. Later we shall explore schemes which promise 


some improvement over the existing methods. First new 
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expansion sets and their convergence properties will 
be examined. Then integral transform techniques will 
be considered. 

Through the treatment of "the particle in a dipole 
field" the unresolved problems to this point will be 
made apparent. 

A final effort to devise an algorithm which would 
eliminate these convergence problems is investigated. 
However at this point the choice of two non-linear 
parameters resulting in maximum convergence remains 
unknown, except for the two limiting cases R->0 and 


R>o where "R" is the internuclear separation. 
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The Schroedinger Equation and Separation 
of Variables 


It can be shown that the Schroedinger equation 


pet ~ > _ > 
{-45V (Z)/r,) (Z,/r5)} Bia) ee ee Jka 1) 
in Hartree's atomic units lk) 


(*) 
can be subjected to separation of variables if use 
is made of confocal elliptic coordinates. 
ibe ry and ro in [{l1] denote the distance of the 
point P(r ,+05r6) to the foci and ¢ denotes the angle 


around the major axis of an ellipse containing P, the 


coordinates of this point can be expressed as 
en =k = 
A s(r,t+r,)/R , w 2(r,-r5)/R [2] 


WoLeEIng tneswave munctiom Wen (jas a product of 


functions each depending only on one individual variable 


WX Mey @ade es LL (Xy) co Mt) ug) (3) 
the DE in [1] becomes 


Le 1 


{ (d/ad) (A2-1) (€/ad) + (d/dyu) (1-7) (d/du)t+[(A27-1) 7+ (1-2) 7] 


(d*/ao*)+R(Z)+Z5)-R(Z,)-Z.)u-p? (A*-u*)F YAU, db) = 0 


where Dp i= See RR: [4] 


(*) Hartree's atomic units are used throughout in this 
Lexi. 
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The last DE can be separated into a system of 


three simultaneous second order lin. (linear) DE's 


[(d2/do*) +m leo) =90 [5] 


[ (d/dd) (X21) (4/Ad) -m? (X?-1) TR (2, +Z,) A-p?A?-A, L(A) = 0 
[6] 
[ (d/dyu) (1-u?) (a/ay) -m? (1-u?) “T-R (21-25) utp? un? +A, 1M (u) = 0 


6 = = A fd 
with three separation constants A,r A, Andam. -orsa 
proper solution we require A, =A5- 


For [5] the solutions are found by inspection 
6(¢) = exp(imd)/V27 [8] 


Since the separation constant m in [5] appears only 
as m* we find that for energy states with mz0 we will 


4 iaene 


have a twofold degeneracy corresponding to m=+|m 
homonucl. (homonuclear) molecule ions where we have 
Z,=Z, gles ihalels Gelchante alg, yl Sey TVA) Aigctatoeleleien 6 
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[ (d/du) (l-w?) (afd )-m? (1-2) © Phone) “eee NG) Te [9] 


This DE is then invariant with respect to an 
operation replacing u by -yp. Hence [9] has a pair of 
lin. indep. (linearly independent) solutions of which one 
is even and one is odd. 

Furthermore the solutions to [9] do not explicitly 
depend on the internuclear distance which we write as R. 

Extensive tables for solutions of [9] have been 


published by Stratton et eh, Se! 
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Solution of the p-Equation 


Restricting ourselves to the treatment of the 


i the 


homonuclear case, according to E. A. Hylleraas 
eigenfunctions of [9] can be expanded in associated 
Legendre polynomials satisfying the following DE : 


cm 


Myer ety ara Roe ety eee [10] 


For the eigenfunctions of [9] we write the expansion, 
M(Zemeu) = E ea) [11] 

The summation in [ll] ranges from m to ~ if M and 
@ are to have the same parity and the summation index 
ranges -i2c0m nt) Sto soa bhempanities of M and’ ¢ diitfer-. 
Furthermore this summation index runs only over even 
integers for bonding states and exclusively over odd 
integers for antibonding states. The symbol 2% in [11] is 
commonly called an "angular quantum number" or "y-quantum 
number" and the expansion in [11] itself in this report is 


generally referred to as "angular expansion". 


Eliminating in [9] the part which is satisfied by the 


associated Legendre polynomials we can write, 
2 es 
ac lis u +A, 2% ( (2+1) Iss (u) = 0 eon 
Q 


For the associated Legendre polynomials the following 


recurrence relation holds, 
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+[{{(2+1-m) (2+1+m)]/[ (2241) (ap penabie (a) 


AUGER) sb YE CL MONE Boe Ty [13] 


Due to the separability of variables the nodes of 
the wave functions have simple geometrical shapes. For 
the homonuclear case the nodes corresponding to the case 
M(u)=0 are two-sheeted hyperboloids with two sheets 
coinciding in the plane perpendicular through the center 
of the internuclear axis for the antibonding states. 

Substituting the recurrence relation for the 
generalized Legendre polynomials shown in [13] into [12] 
the expansion in [11] is defined by the following three- 


term recurrence relation among the expansion coefficients 


[[ (£-1-m) (2-m)]/[ (22-3) (28-1) ]]e,_, 
+ [[(Aj-2 (241) ) 1p 7+ [ (41m) (L414) 1 /[ (2043) (2841) ] 
+[(&—-m) (24m) ] /[ (2241) (22-1) ] Ic, 


+[[(2+2+m) (£4+14+m)] /[ (2245) (2243) ]]e5, 5 = 0 [14] 


This represents an infinite system of lin. homogeneous 
equations which have a non-trivial solution only if its 
determinant vanishes. This determinant does vanish for 
appropriate choices of the separation constant Ans Ze 
wie A. in [14] represents the eigenvalues of the 
corresponding tridiagonal matrix. For Similar expressions 
in the following sections this role of the separation 


constant will be assumed to be stated implicitly writing 


an expression of the type [14]. 
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Since M(u) represents an expansion with infinitely many 
terms a crucial question at this point is whether this 
series converges sufficiently fast for practical purposes. 

Hylleraas, who proposed this expansion, using a power 
series expansion for the separation constant showed that 
for internucl. (internuclear) distances of chemical 
interest the rate of convergence is very rapid. ey 

From the yrecurrences relation [14)7 it vs readily 


apparent that, for small internucl. distances, R=0 


the eigenvalue spectrum of the separation constant is 
A. (2) ae ( 21) Qa enO plas 2 pk enees fabsy 


and since for R=0 we are approaching a hydrogenlike 
united atom we have a rough estimate for the electronic 
energy at hand 
aires 2 2 
E,(R) = (2, +2.) VaK eile) 
Tee 2 oy ee (in Hartree's 


atomic units) [16] 


This limiting behaviour of the electronic energy and 
the separation constant is of interest in respect to the 
application of the computer program and will be referred 


tor Later. 
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Qualitative Inspection of the Rate of 


Convergence for the yw-Expansion 


Considering equation [10] which is the DE satisfied 
by the generalized Legendre polynomials and comparing it 
with [9] where we have in addition the term epee we find 
that the rate of convergence will be slow if this term 
becomes dominant. Recalling that, 

p?= SER Ba? Gye 
the rate of convergence is expected to be slow for large 
internuclear distances. The expansion set constitutes the 
set of eigenfunctions to [9] for the limiting case R=0. 

For R40 the rate of convergence for large m quantum 
numbers is expected to be superior to the rate of convergence 
for small m's. This trend is predicted because for large 
m's? the term “m-/(1-u)" which is satisfied by all terms 
in the expansion, constitutes a more Significant term of 
[9] than if m is small, and it is absent alltogether if m=0. 

Tniisearoumiecimt pal Clough sOt iIstriCchiyeaInturcive NatuLe was 
found to hold when it was subjected to tests in numerical 


calculations. 
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The Radial Expansions 
(A) The Hylleraas Expansion 


We recall the DE for the radial wave functions 


OP Ty 7a) Pana one y (nie /( a=) ) OR pe ALL ( 2) = 6 


Q= R(Z +25) /2 (185) hee) 


il 
The asymptotic form of L(\) for large X} is easily 


found. Dividing (1el by AW 


(ES CURA do /dn (27) (d/d yam Oh Olean 


4(20/X)—p*— (A, 727) JL) = 0 [19] 


Then for large A's the terms with negative powers 


of } can be disregarded. This gives us 
laa /d ji) =p 1 (a 0 [20] 


which iS a second order lin. DE having the char- 


(characteristic) polynomial 


Popo) 0 (21) 
and therefore the homogeneous solution to [20] can 
be written as 


TO) eS exp (-p\)+c,exp(-pr) 2 2) 


il 
If we take p>0O the second term in [22] does not 
contribute to the physically meaningful form of the 


ASYMOUEOCIC ft Orm Of aby oer c5=0. 


(x) This. substitution noldas: throughout and foreune 
symbol eke” used later we have 0,-=Q- 
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When attempting to solve [18] we first investigate 
Suitable transformations to remove bothersome terms. 


The ticst such transtormation is 
Tee ect One [23] 


removing thentecme my (\ —1)" 
One nie ViNe yeh emda dn) — (one 


—2O), FAS.) lai A)  =8 0 [24] 


1 
Changing the variable 


AX = Mees) ' Ae [1,@) ' e107.) (225) 


(24 transtornmss to 
{x (x+4p) (d*/dx*)+2[ (m+1) (x+2p) ] (d/dx)+[- (x/2)? 


+ ((Q/p) -p)x+2Q-p* +m(m+1)-A, ] } Vite = 1 [26] 
For R=0 we have 


ed 4 a2 oly +Zp)e/ (20°) and hence: /i(0/p) =n 27 


uh 
Substituting these limiting values for R=0 into 


[26] one obtains the DE 


{x2 (a2 /dx2)+2 (m+1)x (d/dx) + (- (x/2)tnx 


—~2(2£+1)+m(m+1))} y(x) = 0 [28] 


woiach has a solution 


PH Fe rll 


(Um) 
< n+2 


For ( nie adit) = exp (-x/2)L (x) [29] 


where L is an associated Laguerre polynomial. 
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Hylleraas expands the eigenfunctions of [18] in terms 


of functions corresponding to the limiting case 
x>0 , R*~ and equivalently p-»~, 


defined by the generating function 


pe Re Ginnie) cea ace eee) 
n=m+1 
Ce a / DOr Vat) et) *) [30] 
which ware, solutions of ther DE 
{x (d2/dx?)+(m+1) (d/dx)+ (- (x/4)+n-((m+1)/2))} £(x) = 0 
eo 
Then using the expansion 
y(x) = See cf (n-1,m;x) [323 
andmsubsei tion ett alnto. thesD rea lzZ 6 
or emi iitmelysetd ax) |b (O7/ pn ti (meta) a2) |x 
n=m+] n 
+ [2Q-p*-A, tm (m+1)-4p(n-( (m+1)/2))]} Bin ema ea 0 
(sie) 


we obtain the above first order lin. DE. The second 
derivatives were eliminated making use of the diff. 
(differential) operator satisfied by the expansion 
LuneeiOns. 

It can be shown that for the expansion terms the 


FOoLLoOwing sholds: = 


xf(n,m;x) = (2n-m-1)f(n,m;x)-(n-m)f£ (nt+1,m;x) 


-(n-1)f (n-l1,m;x) [34] 
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1S 
x (d/dx)£(n,m;x) = -—((m+1)/2)£(n,m;x) 


+ ((n-m) /2) £ (n+1,m;x)-((n-1) /2) £ (n-1,m;x) 35) 
Substr tution tor these im [33] leads@to, 


Ee, {£(n-1,m;x) [-(n-1) ((Q/p) -n¢m+1) ] 


n=m+1 
+£(n,m;x) [ (2n-m-1) ((Q/p)-n)+2Q-p*-A, 
+ (n-1) (m+1) -2p(2n-m-1) J 
he(n+l,m:x) [—(nam) ((0/p)—n)]} = 0 [36] 


TMhenvcollecring lermsvwith 1 dentical f(n,m:<) (361s 


obtained in the form, 


Z {[-(n-m-1) ((Q/p) -nt1) ]e,_ + (2n-m=1) ((Q/p) =n) 
n=m+1 
+2Q-p*-A,+ (n-1) (m+1)-2p(2n-m-1)]c 


n 
~[n((Q/p)- (n-m) )]c_ 4 }£ (n,m;x) =O. Loyal 


In order for the expansion in [37] to hold we simply 
demand that all) coerricirents or =the £(n,m;x) ss vanisn. 
This requirement leads to a recurrence relation among the 
expansion coefficients. In our case we arrive at the 


following three-term recurrence relation: 


~{ (n-m-1) ((Q/p) -nt+1) Je, _,+{ (2n-m-1) ((Q/p) -n) +20-p* 


-A,+(n-1) (m+1)-2p(2n-m-1)]¢,-[n((Q/p)-(n-m))]c,,, = 0 


* Hylleraas in his original paper lists this formula 
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Thissexpression sis of; the same typeras, [14] 7and can 
again be interpreted as an infinite algebraic eigenvalue 
problem on grounds of the arguments brought forward there. 

Since the tridiagonal matrix corresponding to [38] is 
not real symmetric (or hermitian for that matter) we are 
not guaranteed that, for the practicable case of truncated 
expansions, the eigenvalue spectrum for the separation 
constants be confined to the domain of real numbers. 

That in the limiting case of infinite expansions the 
separation constants are necessarily real is apparent 
since [6] is self-adjoint. 

Since according to the power series expansion for 
A given by Hylleraas, for finite R the A's should be 
in general non-degenerate, there should be no multiplicity 
of the roots Gb sthevchar. “polynomial of the determinant 
to be expected. Therefore homing in on the roots along 
the real axis should generate approximations to the 


true eigenvalue spectrum rather efficiently. 
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Convergence Properties of the Hylleraas 


Expansion 


Numerically the convergence is found to be best for 
large m quantum numbers, which of course requires high 
& quantum numbers. In addition convergence turns out to 
be better for large internuclear distances than for 
small ones. This behaviour can be expected since the 
basis functions defined by the generating function in [30] 
anevergenfunctions of [26] for the case kes, contrary to 
the angular expansion which is accurate or diagonalized 


for the case R=0. 


(B) The Jaffé Expansion 


If we apply to [6] the transformation, 
HOM (ile exp (pi), (2) [39] 
we obtain the DE, 


{(A2-1) (a?/dA7)+2[ (m+1)\-p(A*-1)] (d/dd) 
+2po\+m(m+1)-p*-A, }L, (A) =a [40] 
where in this section, 

o = (R/p) ((2,+2,)/2)-m-1 [41] 
Jaffé applied to this DE the bilinear or Mébius 
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ee) col) [42] 


which maps the area Re(iX)2>=1 onto the disk with 
radius= tangent to the imaginary axis at the origin 
in the complex plane. By this transformation the area 
around the nuclei is moved to the origin, i.e. the center 
of expansion and the singular point A=(-1,0) is moved 
to the point at infinity in the extended complex plane. 


Thismeransrormse 40). into 


{& (1-—)7D*+[ (m+1) (1-2 )+ (1-&) (m-1)-4p&]D 


+[m(m+1)-p?~A,+2po ((1+E)/(1-E))]} LL (E) = 0 [43] 


The last term in [43],which becomes singular at é=l, 


can be removed by the substitution 
La(é) = 94i(1- eb) 1 CE) [44] 
Ah 2 
WitChetrans corms. 45). nto 
{&— (1-—&)?D7+[? (1-m-20) +2 (o-2p-1)+(m+1) JD 
+o (atm) &+(o+m) (m+1)+2po-p*-A, } Lode ee [45] 
Now if we assume for L, a power series expansion 


bo(E) = 2 aE [46] 


using standard procedures we obtain the following 
three-term recurrence relation among the expansion 


coefficients in’ [46] 


G,-1 | (n-l-9) (n-m-1-9) ]+q, [2n (o-2p-n) + (o+m) (m+1) 
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It was implied by Teac ale that his expansion 
Ps 2 (m/2) O ne n 
UO (ON ah) (A+1) SPS qi ((A=1) / As1)) 
[48] 


converges faster than the one found previously by 


BoA. plies 


However, we can show that the 
eigenvalues of finite matrices of equal size generated 
by truncebionpor the infinite set of lin. (linear) 
equations of either [38] or [47] are identical. 

In other words, considering an equal number of terms 
in [38] or [47] generates identical finite eigenvalue 
spectra in both cases. Having an equal number of terms 
generates in both cases a quantitatively equal degree 
of approximation to the exact infinite eigenvalue 


spectrum of the separation constant of the radial 


expansion. 


# Proof : 


We wish to evaluate the roots of the char. 
polynomial corresponding to the finite algebraic 
eigenvalue, problemidet (B-Al)—0.) Injoun case 6 is a 


tridiagonal n-dimensional matrix with 


2 : leny oe = : 
Toad Spee ee cise! een y 


[49] 
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The elements of Bewould be® the@coctiicilentsminmc ne 
systems of lin. equations of [38] or [47] respectively, 
excluding the separation constant A in the diag. (diagonal) 


terms] Then therchar? polynomial 
De A)eeedecE EAT) [50] 


where I is the n-dimensional unit matrix, may be 


evaluated as follows 


Let Py; (A) denote the leading principal minors 
of WBeAL)sotsorder k. Then expanding Py (A) by its 


last row, we have 


Watchin Ke = 23 7 Aip staterete 1) 


and antvtially, Po (A)=1 ; P, (A) =a,-A ay 


This means that the char. polynomial Py (A) 
depends only on the diag. (diagonal) elements and 
the binary products of symmetrically positioned 
pairs of off-diag. (off-diagonal) elements. 

One can readily see, that in [38] and [47] the 
diag. elements are equal, if it is kept in mind 
that the expansion index in the Hylleraas expansion 
begins with n=m+l whereas for the Jaffé expansion it 
starts at n=0. In addition, keeping the difference 
of the “expansion indices ain" (S8)]and [47] in mind, 
it can readily be shown that the binary products 


of off-diag. elements in [51] are equal using either 
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Hence according to what was said previously 
about the evaluation of the char. polynomial of 
tridiag.(terai.agonal) ’matrices,. the char. 
polynomials have to be identical in both cases 
and of course they will have the same set of roots. 
This means that the finite eigenvalue spectra for 
the separation constant are identical in the case 
of either a truncated Hylleraas or Jaffé expansion 
as long as an equal number of terms is considered 
in both cases. 

Bree end of proof # (x) 

The Jaffé expansion was used in almost all the many 
calculations of this type to date. However it was just 
shown that it provides no improvement in respect to 
convergence properties over the previous Hylleraas 
expansion. The major advantage of the Jaffé expansion 
probably lies in the fact that the physically meaningful 
domain of the expansion is a finite interval, making it 
easier to handle when the wave function is to be used to 
extract quantum mechanical properties of the system using 


numerical integration techniques. 


(*) Using Laplace transform techniques it can be shown 
that a functional relationship exists between the 


Jaffé and the Hylleraas expansion. 
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CHAPTER I] 


In this chapter expansions not previously known 
are presented. It is believed that the one-particle 
two-center problem is most efficiently treated making 
use, or thee theory of special functions. The theory of 
the generalized hypergeom. (hypergeometric) series 
and its degenerate cases is used extensively in this 
chapter. 

We often do not know in advance whether the radial 
or the angular solutions of the GSW-DE are generated 
by a proposed expansion. Often expressions like 
"radial expansion" or "angular eigenvalue spectrum" 
etc. are used in this text. These refer to the type 
of eigenfunctions and eigenvalue spectra generated by 
the corresponding expansion. For an expansion to be 
an angular expansion the wave function must remain 
bounded in the range |u|<1 and for the radial 
expansions the wave functions have to remain bounded 
for A21l. For the angular expansions the homonucl. 


(homonuclear) case is treated separately. 
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Expansion of the Angular Generalized Spheroidal 


Wave Functions in Ultraspherical Polynomials 


The angular GSW-functions are power series 


expandable solutions of the p-equation. If we use a 


power series expansion for these functions the 


homonuclear case can be treated separately. 
If we apply in [7] the transform 
We 
M(u) = (1-u2) 7/2, (u) 


the following DE is obtained 


{ (u?-1) (d?/du*)+2 (m+1) py (d/du) -p*u?+Ru (2, - 


+m (m+1)-A, } M, (u) = (0 
ae Z,=25 a power series expansion of M, 
Ze n 
M, (hu) = a qh 


[52] 


>) 
[3] 


[54] 


generates the following three-term recurrence 


relation among its expansion coefficients 
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G49 { (ntl) (n+2)1 = 0 


0, as Ger ateraiers for bonding states 
(gees ARES as wt PR PURI fonrantsbondangesstates 1) (55)] 


For Z,+Z. thes transformation 


M, (u) = exp(-py)M, (i) [56] 
yields the DE 


iG 1s) De 2.1 (mt) pa DER ZZ 2p (mal) 1 


+m (m+1)-p*-A,} My (1) ae) (57) 


A Taylor series expansion for Mo 


M,(u) = B salen) [58] 


is defined by the three-term recurrence relation 


among its expansion coefficients : 


nop b2 (ntl) (ntmt1) ]4+q, [n (n-1)+2n (mt+1-2p)+2po+m (m+1) 
Baad ay WP (nH l oer mO IIR (2 7e-45.))/e(2p).1 (ms) 


[59] 


For the homonuclear case this expansion can be used 
as well. However we are going to show that rapidly 
convergent orthogonal polynomial expansions can be found 
for the angular eigenfunctions leaving little practical 
importance to the power series expansions in [54] and 
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A. The Homonuclear Case 


The expansion in associated Legendre polynomials 


found by Hylleraas ‘?) 


is very rapidly convergent. 

For an expansion in ultraspherical or Gegenbauer 
polynomials this is the case as well. How the Legendre 
and the Gegenbauer expansions compare in their rate of 
convergence for different states and internucl. 
separation was never systematically investigated by 
the author.The reason why a new expansion for the 
homonuclear case will be presented here is because 
uSing this type of expansion, we find that we can 
proceed rather easily to the treatment of the 
heteronuclear case, for which no expansion has ever 
been published with a comparable rate of convergence 
and ease of handling. 

The GSW-function to be expanded for the homonucl. 
case is M, satisfying the DE [53]. 


We try an expansion, 


= (a) 
M, (u) = Pedy CG (uu) 


(a) 


where the Cy are ultraspherical or Gegenbauer 


Dolynomialsee indicating sthe ditt.) Operator satisiied 
by these we get our notation straight 


Greeti = { (u?-1)D7+(2o0+1)D-2 (2+2a) } eae = 
D=d/du 


ee 
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[60] 


0 
[61] 
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[61] is a hypergeom. DE with three regular 
Singularities. 

In the process of setting up the expansion we will 
make a choice for the parameter oa such that a bothersome 


term will disappear. 


The ultraspherical polynomials satisfy the recurrence 


relation, 


ute (1) = (241) C09) (u)+(2t2a-1)c'%) (u) lV 
[2 (+a) ] ee 


from which we derive the desired property, 


ec ier) = {£ (241) (242) /(Lta41) 105%) Cu) 
+EL (241) (2420) /(Ltat1) 142 (2420-1) /(t4a-1) 110)" (y) 


#£(2420-1) (242 (0-21))/ (240-1) 166%) (Wy }/L4ta)] [63] 


Now we set up the expansion, 
(u) 


Sq, {[2m-2o+1] wD-p?y2+m(m+1)+2(L+2a)-A,} i 
[64] 


= 0 
where use has been made of the diff. operator Go 
satisfied by the ultraspherical polynomials. If for the 
parameter a we make the substitution, 
a = m+(1/2) 
PniewcoctELcouvent. tony tresturst derivative. ins (64) 
(63) 2ne164)] one 


vanishes. Making use of the property 


gets an expression : 
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; Gy El-p? (241) (£42) 1/[4 (2a) (24041) 11099) (u) 


Hier eA) i Cee + 2a) eros 
+[ [2 (420-1) ]/[2+a-1] ]] +m (m+1)+2(2+20)-a,}ch™) (u) 
+[[-p? (2420-1) (242 (a-1))1/[4 (2ta) (2+a-1) 11 Ch) (uy } 


= 0 [66] 


Collecting equal polynomials we get from this a 
three-term recurrence relation for the expansion 


coerfiicirents: in [60], 


dp_ol [=p (241) 1/14 (24+a-2) (240-1) 1] 
tq, [[-p?/(4(2+0)) ] [1 [ (2+1) (£420) ]/[2+0+1]] 
+[ (2 (2420-1) J/[<t+a-1] ]]+m(m+1) +2 (2420) -A,] 


+q [[-p? (2+20) (2+20+1)]/[4 (2+a+2) (L+at1)1] = 0 
Q+2 167] 


with re a rane iy ate OPC for bonding states 
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B. The Heteronuclear case 


In the heteronuclear case M which satisfies the 


Dak 
DE [57], can be expanded in ultraspherical polynomials. 
In order to shorten the expressions the following 


substitution will be made 
B= R((Z,-25)/2)-p(m+1) [68] 


In addition to the recurrence relation [62] the 


following differentiation formula is needed : 


(a) 


(u2-1) (A/au) Co) (a) = {2 (2F1)CF% Cu) 
~ (2420) (2420-1) C)%} (u) 1/12 (ata) } [69] 


Prom f5/ 1), making use of the diffi. operator Go satisfied 
by the ultraspherical polynomials and assuming for M, an 
expansion in ultraspherical polynomials, 


(a) 


Mo (u) = A G,Co tu) [70] 


one gets an expression, 


z q, {[ (2m-2a+1)u-2p(u*-1) ] (d/du) +2Bu 


+2 (2+2a) +m (m+1) —p?-A, } Geo ee [71] 


The parameter a is chosen such that the linear term 


for the first derivative vanishes, 


P= 2 Ot ie eee nt fl / 2.) L721 


which is the same as used for the homonucl. case. 
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Substitutions .or thewproperties®. [62] .and (69) sorrtne 
ultraspherical polynomials gives an expression of the 
same type as [66]. If we collect equal powers in this 
expression and require their coefficients to vanish we 
arrive again at a recurrence relation among the expansion 
coefficients in [70]. It is the following three-term 


recurrence relation 


[[2(B-(2-1)p)]/[2+a-1l]q,_4 
+[&(24+2a) +m(m+1)—-p*-A,] aq, 


+EL(2Q42a) (Bt (£42041) p)]/[&tatl]]q,.4 at £8 [73] 


where = Oy gleg 2 peciperr ene es ® 


Examining this recurrence relation we find that for 
the expansion [70] the corresponding algebraic eigenvalue 
problem is diagonalized if R=0 and p=0 i.e. in the 
BUnLeeGasatom. slim. 

BLOM LLStingseiim Piguet. 2t scan hes seen  ehatetaking 
only 15 expansion terms a reasonable number of significant 
figures is obtained with this exvansion. The eigenvalue 
referred to in Fig. I is the second lowest eigenvalue of 
the angular expansions which is listed in Appendix II. 

For Hell the state Inet,m>=|6725 > was calculated 
for five hundred internuclear distances from R=0.2 to 
R=10020 with tansincrement of AR=0. 2) to: tls “significant 
figures using 40 expansion terms for both the angular 


and the radial expansion. 


: % oe & iw 


; ‘ oe Ty 7 7 
7 me . ’ ‘ ‘ i ow i ; nae PS 
- u | ; is #° © iit ae faa 
me ai es i 


eas phate hat 


7: Yaad 
Wit * ei f iain i lo Phes\ we) a a aes 


ie P 
€ ¢ Saceey >! = 1%, Ay J BoC) 5 £ ties sd bes ee Pa a ; a 
_ ¢ 
snlun ade Ut® Grogs) Cos teres (nei tues na FO lel 


= * ela A Fal > 
j {-s waite lio | al a) ea. . | e a | * : lj al 


4 ies, @ 


7 . 
h riya I -ofn kere apets! 
fy da 
i? gh 4 J a 
« 5 a 
‘ ‘ 
7 iow 
’ Tat YH8.-¢ ty yt s i ntase 
cine + ae Like bee ees {AR 
aid ay t Spal Gat? SToe : = werreage | ee 
' of 
ana stns 4449 Ata esd pas AT V+ ge" 
She | 
‘s ION4 arny Ps ry 5 D Ul Us & 7371594°.% - J "Ce ; [GAM « 
\eiclavasefe «ay ..09 Aeeaae .t ied | niet (tte st 
u | " : 7 : 
_ a 7 ioe” =e on a 
i] if Ve <“sf1 Cine. 8 mes: ait ee hi has ’ N, = 9s oo ~? 


f . ae sd <xfovah AL rset =). tee hae Paste | beim eas sim as 
a - 
Poe: 4" “plates Br «1 \s eel ‘adegins 
on 

ae 4. ver in Tee ail 


ia Te ow b> 
i Precis ae aire bat 


= 
- 


: 
iad 


9 


28 


The computer program used for this purpose is 


listed in Appendix I. 


Uke ee Bees The number of significant decimal 


(*) 


figures for the second lowest 
eigenvalue (Separation constant !) 
with m=0 of the angular expansion in ultraspherical 


polynomials using 10(Fig. Ia) and 15(Fig. Ib) expansion 


terms. 

(a) (x*) 
jake Ome 2 3 4 5 6: cal 8 9 10 

Q, 
0n5 Som 2 i Se Gy a he 4 3 
eS PL Om 2 7 6) 54, 4 3 3 3 
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9.5 hee ae wa ayes) OS eA. aks 


(*) Double precision was used within a Fortran IV source, 


which means a mantissa of 14 hexadecimal figures. 


(**) Qj=4R(Z-Z>) 
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For the angular expansions one can take, 


(i) i ie a 
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which essentially corresponds to a change in sign for 
the variable and replacing p by -p in the transformation 
[56]. For the expansions in ultraspherical polynomials 
the rates of convergence in respect to the separation 


constant are equal in all four cases. 
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Expansion of the Angular Spheroidal Wave Functions 
in Spherical Bessel Functions 
It can be shown that for the angular spheroidal 


wave functions (homonuclear case) an expansion in 


modified Bessel functions is possible, 


M, (11) qu, (Pu) q,, (pu) "K, (pu) [74] 


2) 4) 


M, (u) qv, (Pu) vr (pu) K, (pu) [75] 


53M 8M 
H] 
5M 


where the K's are modified Bessel functions of the 


third kind and M, is scetined by [52 and [53] 


The linear differential operators satisfied by the 


expansion functions are, 


{xD*+ (1-2n) D-x} u(x) 


{xD2+(14+2n)D—x} vy, (*) =7/0 (a7 


Il 
Oo 


[76] 


Examining these two operators we can easily see 
that the u's and as are either even or odd functions. 
The properties needed for the expansion are easily 


derived from the theory of Bessel functions, 


xDu,, (x) = énu, (x) -u i.) Oe) [78] 
D*u, (x) = (1-2n)u__, (x) tu, (x) [79] 
xDv (x) = Sper ot) ravines 30) [80] 
D*v_ (x) = (1+2n)v_ 4 (x) tv, (x) [81] 
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Scaling the variable, 


[82] 


* 
Il 
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[53] becomes, 
{ (x*-p*)D*+2(m+1) xD-x*+m (m+1)-A,} ear [83] 


Substituting the expansions [74] and [75] into [83] 
and making use of the linear differential operators [76] 
and [77] satisfied by the expansion functions one obtains 


the following expressions : 


y q_{-p?D?7+[2(m+n)+1]xD+m(m+1)-A.} u_(x) = 0 [84] 
ret 2 n 
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5 r,{-p*D*+[2(m-n) +1] xD+m (m+1) -A, } Vv, (x) 


The ditrerencial relations: (801 and -[S1]"for Che 


Yes 


nS can be obtained from the ones for the UL 
replacing n by -n and the same holds for the differential 
equations for the u's and vss Therefore the recurrence 
relations for the two different expansions will predictably 
differ again just by the sign of n. 
Hence we stop working on the expansion [75]. 
Substitution of the differential relations [78] and 
[79] into [84] and applying standard expansion techniques 


gives rusestosthne.tollowing, three-term recurrence 


relation among the expansion coefficients in [74] : 


q-1 [1-2 (ntm) ] +a, { (2ntm) (2n+m+1) -p *-Ay] 


Sealey (URES [86] 
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The algebraic eigenvalue problem associated with 
{[86]-is bisected at n=+. 


If we take 


or Tle Server eterateme hr Kany jets etere ef) det eres e (ea) ke sie 
-o<k<o [88] 
the eigenvalue spectra for the homonucl. angular 
expansions are obtained. The expansion with the 
expansion index as in [87] is more rapidly convergent 
than the one with the expansion index being an integer 
number as in [88]. 
The Bessel functions of integer order can be used 
for even or odd angular states, whereas for the 
spherical Bessel functions the following holds : 


TL SiO ieee cipanaye even angular states 


TL en vara odd angular states 
If we seek an expansion in modified Bessel functions 


of the first (and second) kind 
Bs n 
M, (u) = z q,, (pu) I, (pu) [89] 


the following recurrence relation holds for the 


q, sin Pres) er 


-p? [2nt+1] q,,,+[(2n+mt1) (2ntm) -p*-A,]q, 


+[2(ntm)-1]q,_, = 0 [90] 


If we choose expansion terms where a negative power 


of "pi" was factored out the resulting change in the 
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recurrence relation for the corresponding expansion 
cOerricientsawill only be a change: in signyroren: 
Considering [90] as a system of linear equations 
it would differ from the one in [86] only by the sign 
of the off-diagonal elements. This means that for an 
equal number of terms considered their eigenvalue 
spectra will be identical. This can be shown to hold 
making use of the same arguments used to exhibit the 
equivalence of the Hylleraas and the Jattc expansion. 
The equivalence of [90] and [86] in terms of convergence 
properties of Ay however is even more obvious, because 
we can show, that a pair of equidimensional matrices 


whose off-diagonal elements have opposite signs but are 


equal in magnitude, have identical eigenvalues. 


Theorem : The two equidimensional matrices B and B® 
differing from each other only by the signs 
of all their off-diagonal elements (but not the diagonal 


elements !) have identical eigenvalue spectra. 


COO ieee Imagine the determinants of "det (B-1I)=0" 
and “det (B°-eI)=0" factored out inn! 
different permutations of the matrix elements, 
appropriate to calculate the determinants. 
Since every one of these permutations contains 


n-l off-diagonal elements the resulting 


characteristic polynomials will differ only 
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(n-1) of ally their 


byl aslac tormof 7 (-—1)) 
coefficients. Therefore they will have the 
same set of roots. 


eng ot /proon 1 


For R=0 the eigenvalue problems associated with the 
modified Bessel function expansions are diagonalized. For 


example taking [74] the analytic expression for R=0 would 


n= | 2-m| k=n 
= ees i 7 [[A,- (k-m) (k-m+1)]/ 


pko2 (malo) b(pnlien aoe 


n= (1/2) (PH) } 
(324 
where aS mentioned previously the separation constants 


for R=0 are : A,=h (X41) eB 
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Expansions of the Generalized Spherical Wave Functions 


in Confluent Hypergeometric Functions 


When in the following two sections we refer to the 
GSW-DE we usually assume it to be in a form where the 


coefficients are specified as follows : 


{[a,xta sa) Dee Lb b x+box°]D+[k +k, x] } Sw(x) = 0 [92] 


2 1 


and in general we assume that in [92], 
lel 4 Sonia [93] 


We are going to show that for any choice of the 
parameters in [92] an expansion in terms of any solution 


to the linear differential operator, 


Ly 9, (OrYFX) = {xD*+ [y-x]D-a}g, (a, ¥;x) (9) [94] 


can be found. 
The power series expansion of the above DE is 


Gy (aryex) = (PCy)/T(@)] EF  CEP (atm) /T (ytn)] [x/n17} 
n= 


OU, ¥ox) [95] 


Using the notation of the generalized hypergeometric 
series, 
J, (a,yixX)=,Fy (a,yix) [96] 
which we denote by Humbert's symbol ¢(0a,yY;x). 
Theerocopsecmepne mnalCcid meequamuon Of [94] eares0,and 
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TLE we pul in. [94] 


g(arvixd=x tn (oryex) are 
we obtain the DE, 
{xD*+ (2-y-x) D+ (y-a-1) } hy (a,¥:x) = 0 [98] 


which is of the same form as [94] and hence its power 


series expansion according to the adopted formalism is 
hy (a, yi) = $(a-y+1,2-y;x) (99) 


Therefore (for non-integer y's) a second linearly 


independent solution is 


gylaryex) = x97Y o(a-yt1, 2-¥:x) [100] 


The functions Dy and 35 form a fundamental set of 
solutions of [94], provided that y is not an integer. 


SUDSCI Uinta, 94), 
Gy (a,yix) = exp(x)h, (a,yix) [101] 
fi eCransc OOMS, LiNbO: 
{xD’+ (y+x) D+ (y-a) } hy (a,Y7x) = 1@ [102] 
Replacing =< by s-xag1Ves, 
{xD*+ (y-x) D+ (a-y) } shtml = [103] 


which again is of the same type as [94] having the 


following system of linearly independent solutions 


®(y-a,¥Fx) ae 


hy (ary? -x) 


h,(a,y:-x) = iy Dino ete hee) [105] 
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TMereLore. according to tnersupstitution | LOL], the 


solutions to [94] can as’ well be written as 


G3 (a,¥;x) Exp (x) O(V=O7 Yk) [106] 


Gy (a, Vix) Seniesa eer 8) [107] 


Since [94] can have a fundamental set of only two 
lin. indep. solutions and considering the limiting 
behaviour of S14 at the origin, one comes readily 
to the conclusion that J) =93 and 35=Sy- 


The relation 


OA Oray ae) Ga EXD (X) Oyen 73x) [108] 
is indeed known as Kummer's transformation or Kummer's 
first formula, which in turn constitutes a limiting case 
of Euler's transformation of the hypergeom. series. 

Frequently the solutions to [94] are expressed in 
terms of ¥ functions. The ¥ function is readily obtained 
by seeking an integral representation of a solution to 
[94] using an appropriate kernel 


Gg(ar,yix) = £1 [exp(-xt)S(t)]dt [109] 
0 


where Zo and Z, are complex numbers to bes property chosen. 


Using standard operations one finds that the kernel S 


must satisfy 
{t({1-t) (d/dt)+[ (2-7) ttl-a] } S(t) = 0 (i104 


and Zo and Z4 must be chosen such that 
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Z 
-exp (-xt)t (1+t)S(t)| 0+ = 0 [aaeTay 
0 
sinees (LUO} isa 'ftirst ordey lin DE at May readily 
be integrated. 
Sey eo are) oe Gaby) 
We shall assume that Re(z)>0 and Re(a)>0O. Then [111] 
will be satisfied for ZQ=0 and Zr Therefore an integral 
representation of the type [109] can be written as : 
Go (a, y7x) = Cee) 


creole) 


= J, Lexp(-xt)t Gee tana [113] 


We can rotate the path of integration in the complex 
plane, i.e. Z4=© is not unique. One can write more 
generally Z,=|~|exp (id), which extends the domain of 
definition = 

wla,ysx) = (P(aj] > gl°le#P 49) pep (ext ye OD) (ret) VOD at 


with Re(a)>0,-1<o<T and -k1<d+arg (x) <7. eke a 


OLSCougse, according (LO stie trans! Ormatlonse!97,|)eand 
[101] to which the DE of the confluent hypergeom. 
functions was subjected, one finds a set Seg 


corresponding to Ji 4 found previously 


Go(a,yix) = ¥ (a, y7x) PLES 
Ge la,yrx) = x TOY y (gay), 2-7) [116] 
Go(a,y?x) = exp(x)¥ (y-a, 7 -x) shy Al 


gg (arvix) = exp(x)x7Y) ¥(1-0,2-y7-x) [118] 
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Atethis point welcertainly would likestosknow first 
of all which ones of Je_g are identical, pairs and in 
addition we wish to find out how the ¥ function 
relates (tO rtheon @eninetion.= 

We start with the second part of our quest stating 


the important relationship, 


¥(a,¥ix) = [T (1-7) /T (a-yt+1) 16 (a, 3x) 
SILT Cys o (eel ao es [119] 


(6) 


which we borrow from literature. 


A readily apparent consequence of [119] is, 


Viasyex) = x VY) v(a-y41, 2-¥7x) [120] 


1 .ee Ge=Je- Then S45 and Tg must be a second linearly 
dependent pair. We find : Jg=—exp (sgn (Im(x) )lym)g_. 

Finally we list all the relationships required for an 
expansion of Sw in terms of confluent hypergeometric (x) 
FuncwLons:. 

Eight solutions to the confluent hypergeometric 
differential operator were presented. Since we can 
choose either aor y as the expansion index sixteen 


different expansions of the Sw's can potentially be 


generated 
Sw, (x) = Y Ty on? (rey F%) = y Senos eos 
fele2 al 
SW) el SIM og Weeue np.) y Xo) Tears) 
2 ner 2 Te ite 2 [122] 
(x) The Sw's are the eigenfunctions of the GSW-DE 
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SWy (x) 


Sw 


10 *) 


Sweex) 


ae 


Sw, 5 (x) 


Sw 


13 °*) 


Sw, (x) 


14 


Sw). (x) 


Sw 


16 ‘*? 


Si) sti ei) eli) eis je) eat sii 


d 
Pe 0 


| 


53M 


shim: tsi Rey ler te) felis “Sri je)i) is) i) yeti 


> 
nN n 43, 


Tan blege nex) i q. 


x6!) vin, ce rx) 


x (nol) 


@(c5,n;x) 


x (not) 


= ¥(Cg,n;x) 


ee pecs x) 


aD, 
e ¥(N,C) 97%) 


P< 
Shige eee 


ax 
LG seen nae 


To Beles one 


413,n 


(nc 


Sly e. aithia 137%) 


(cj 471) -X,, 
aan ey eu(n Gc 


Shun ehiny Curses 


(n-1) _-x 
Ay 5ans e O(c) ,n 


Shins, clits Suge 


(n-1) _-x 
Cie are ee 


eoWe scene eee 


40 


nJ 3 


nda (qr nix) [124] 


y g Gg {nc es) 
n Shawayer is) 5 [125] 


29 Crear, Cras 
Gena 6 6 [126] 


: W797 (C7 enix) [r274 


a Gg nIg (Cgrni x) [128] 


dg nIq(MrCgiX) [129] 
i 910,n910 (B13 97*) 

[130] 
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Now we list the 


xD*+(c,-x)D-n 

xD*+ (c,-x)D-=n 
xD*+(n-x)D-c, 
xD*+(n-x)D-c, 

xD*+ (2-c,-x) D+ (c,-n-1) 
xD*+ (2-c ,--x) D+ (ce-n-1) 
xD*+ (2-n-x)D+(n-c.-1) 
ae (2on=x) Dac ee!) 
xD*+ (Cg+x) D+ (Cg-n) 
Dons, gtx) Dic, Nn) 


) 


2 — 
xD “+ (n+x)D+(n Cia 


) 


2 — 
xD“+ (n+x) D+(n Cio 


xD-+(2=c 3t*) D+ (1-n) 


L 


xD*+ (2-c 4tX)D+ (1-n) 


il 


xD?7+ (2-n+x) D+ (1-c, -) 


hs) 


xD*4+/(2—n+x)D+(1L—-c. -) 
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three-term recurrence relations for 
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[140] 


[141] 
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[143] 


[144] 
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the different expansion sets : 


(c,-n)g, (n-1,c);x)+ (2n-c)+x)g, (n,c) 7x) 


~ng, (nt1,c);x) = 0 


J (n-1,cy;x) -(2n-c5+x)g, (n,cy;x) 


+n (n-cj+1)g, (nt1,c5;x) = 0 
n(n-1)g3(c3,n-1;x)-n(n-l+tx)g,(c,,n;x) 


+ (n-c 3) xg, (c3,n+1;x) = 0 


(n-c,-1)g, (c,,n-1;x) -(n-1+x)g, (c,/n;x) 


+xg,(c,,n+1;x) = 0 


(ce-n)g_ (n-1,c,e;x)+ (2n-ce tx) g, (Nn, Ce; X) 


“ng, (ntl,c,;x) = 0 


g_(n-l,cg;x)-(2n-cetx) gg, (n, ce; x) 


tn(n-cetl)g, (ntl,ce;x) = 0 


n(n-1)xg7(c pn-1;x)-n(n-1+x)g_ (c,n;x) 


" 
+ (n-c1)g7(c7,nt1;x) = 0 
(n-Cg-1)xgg (cg,n-1;x) -(n-1+x) gg (cg,n;x) 


+99 (Cg,n+1;x) = 0 


vc —N) Fg (n-1,Cg7x) + (2n-CgtX) Fg (N, Cg? X) 


9 
“Ng (n+1,Cg7x) = 0 


319g (n-1,C, 97x) ot (2n-C, 9+) dy Q (nN, Cy 97%) 


tn (n-Cy gtl) 949 (ntl, cy gi) = 0 


n(n-1)g,, (Cry piss) ~n(n-1+x)gj4 (Cc), nix) 


+ (m-c) 4) *9}q (Cy, Mtl x) = 0 
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[154] 
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Eksyoe 


TAP eG 


[158] 


[594 


[160] 


LoL] 
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43 
(n-Cy 5-1) 945 (Cy 9en-1;x) - (n-1+x) gj 5 (Cy 9/N;x) 


+XGy 5 (Cy 5,nt1;x) =a [165] 


(C7 3-N)g,3(n-1,c, 37x) + (2n-c, 4+x)g)3(n,Cc) 37x) 


“ng, 3 (nt+1,c, 37x) = 0 [166] 


Sy 4 (n-1,c) 47x) -(2n-cy 4+x) G44 (nC, 7X) 


tn(n-c, 4t1)g, 4 (nt1,c, 47) = 0 i671 


n(n~1)xg,._(c),,n-1;x) -n(n-1+x)g,_(cy5 nix) 


+(n-c (c)o/ntl;x) = 0 [168] 


iis Sais 
(n=c, 5-1) xgy 6 (cy ¢,n-lix) - (n-1tx)g, ¢ (cj) 6 nix) 


tg, 6 (Cy¢ntlix) = 0 [169] 


These three-term recurrence relations can be derived 
from Gauss' relations between contiguous functions. 

The differential relations listed below are readily 
derived from the theory of Gauss' hypergeometric series 
as well. Having these we will have everything required 


for the expansions at our disposal. 


xDg, (n,c, 7x) = n{g, (nt1,c,:x)-g, (n,cy7x)] [170] 


xDg. (n,Cy7x) = n[(n-c,+1)g, (nt1,cy;x)-g (n,Co7x)] 


(n-c,t+x)g,(n,C57xX)-go (n-1,cy;Xx) Ryda | 


Dg, (c3,n;x) = [ (c/n) -1]g, (c3,n+1;x) +g, (c3,n;x) 


(1-n) [9 (c3,n7x)-g3 (c3,-n-1;x) ]/x iy 


Dg, (c,.n;x) = -G4 (Cy ntl;x)+gy (C4 /n7x) [173] 


ef, Jo (n,Ceex) = ng, (nt1,c,.7x)+(co-n-1)g, (n,ce;x) (74) 
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xDg¢ (n, Ce; x) = n(n-cetl)g¢ (nt1,ce;x)+(ce-n-1) g¢ (n,c¢; x) 


(n-1+x)g¢ (n,cgix) -g¢ (n-1,C¢;x) a BG fa 


Dg4(c7,n;x) = [{ (cj/n)-1]} 95 (cz ,n+1;x)+(n-1+x)g7 (cz ,n;x) ]/x 


= (n-1)g_(c5,n-1;x) [176] 
XxDGQ (Cg nN; x) = “Gg (Cg ,ntl;x)+ (n-1+x) gg (Cg,n;x) Faker 
XDJg (N, Cg; x) = NGg (Nt+1,Cg;xX)- (Ntx) gg (N,Cgix) P1733) 


xDgj 9 (NC) ix) = n(n-C, 9t1)gj 9 (nt1,c) 9x) - 
(nt+x) 44 (n,C) 97x) = (N=Cy 9) Gy 9 (Nr C] 7X) 


Clie oc) [2794 


“S10 10 


Dg, 3 (C44 nix) = [(c,,/n)-119,, (¢,, ntlix) 


[(l-n-x)g,4 (Cy), /nix)+(n-1)g), (Cy, -n-1ix)] 
/X [180] 
Dg) 4 (Cj, 5/NFx) = “Jz 9 (Cy ont; x) [181] 
xDg} 3 (n,C) 37%) = Ng, 3 (nt1,c) 37x) +(c)3-n-1-x)g, 3 (n,c, 37x) 


[182] 
XDGy 4 (MCygex) = n(n-cy 4t1)gy, (ntl, cy 47x) + (cy y-n-1-x) 


Dy 4 (MC, gi) = (n-1) gj 4 (M,C) 47%) -Gy 4 (n-1,c) gi) (i334 


Dg, 5 (Cy 5/nix) = [[(c,./n)-llgj.(cyo-ntl;x) 


+(n-1)g,.5(c ,5.nix) 1/x = “Gy 5 (Cy5n7x)+(n-1)g,. (Cc) 5,n-1;x) 
[184] 


xDGy ¢ (Cy gN7X) = “Fug (yg ntlix)+ (n-1)g)¢ (cy ¢.nzx) [235] 


Instead of deriving explicitly the expansions with all 
the steps involved we restrict ourselves to listing the 


essential information about these expansions only. 
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The procedures are formally the same as the ones 
shown in detail deriving the expansions in ultraspherical 
polynomials and others. In the present case a detailed 
listing of intermediate results would generate many pages 
of tedious formulas however. 

We can divide the generalized spheroidal DE through 
by a constant factor such that in [92] we have aj=1 and 
we can scale x without loss of generality such that bj=fl, 


with the sign to be chosen depending on the type of 


expansion employed, 


S Sw(x) = [(a,x+x*) (d*/dx*)+(b)+b,x+x*) (d/dx) 


i 
+(k)+k, x) | Sw(x) = 0. [186] 


Applying to this DE the properties of the basis functions 
of the expansions and furthermore choosing the parameters 


c, such that an undesirable term vanishes, we find that 


k 


we end up in all sixteen cases with a three-term recurrence 


relation among the expansion coefficients of [121]-[136], 
Cc — 
Seow, (2c) anes 2 coe a, (SkIx) = 10) [187] 
lt oa wi) (n) + W2(n) g aye: Ga) 
k sional k k,n k k,n+l 
= 0 [188] 


First we list the parameters Cy of the expansions 


=Cg=Ci9 = bo/ay [189] 


= 2-(b,/a,) [190] 
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[194] 


spheroidal wave functions to which the above relations 


for the expansion sets were applied are the following : 
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The coefficients in the three-term recurrence relations 


as formulated in [188] are 
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Expansion of the Original GSW-DE in Confluent 


Hypergeometric Functions 


If the original GSW-DE is expanded using the 
expansions proposed in [121]-[136] the following 


substitutions have to be made : 


avert 2021 
by = 4p(m+1) 213] 
by = 2(m+1-2p) [214] 
ky = 2Btm(m+1)-p°-A [215] 
kK) = B/p [216] 
where 
B = ae -p (m+1) 
Q. = *R(Z,+Z,) 4 Q, = 4R(2Z,-Z,) [ee 


Independent of the substitutions to be made for 
a rbo Geter, use’ of ‘either Orormy functions: will gave 
rise to equal rates of convergence of the separation 
constant, i.e. identical truncated eigenvalue spectra. 
This becomes easily apparent from the fact that the 
diag. elements of the associated eigenvalue problems 
are equal sin both cases. Of .counse, thesactual wave 


funccions, cdatrer sf) aptunitesnumber,oL, terms, 1s 


considered. 
The constants C116 vat (UEPM Pea RG Way iN Relelelenctelstiates Sete) 
[189]-[194] are : 
CAC =P Cc). SC = mtl [218] 
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For themsiqns toftprandiR ‘themfourscombinations 
shown in [73.5] are possible. Furthermore the 
algebraic eigenvalue problems associated with Ji -J16 
are intersected in two places if m+0 leaving a 
sub-block of size mxm and two additional sub-blocks 
with infinite size. If m=0 only the latter are 
present. Performing numerical work it was found that 


in cases where the expansion index was the second 


parameter of the confluent hypergeom. functions, i.e. 


Sw) (x) = ' CRU celles a ,9) [224] 
or Sw) (x) = z Gao ot AC nie) 22) 


the angular spectra for the separation constant 
were generated. Exploring all of the eight such 
expansions and considering the four sign combinations 
for R and p eight different eigenvalue spectra were 
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with elghte different rates of convergence: An interesting 


fact is that all of the eight eigenvalue spectra can be 


found by taking one individual expansion, considering 


all four sign combinations and the upper and lower 
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infinite sub-blocks separately. For the case m=0 
only four different eigenvalue spectra occur. 

We recall that for the angular expansions in 
ultraspherical polynomials the choice of the signs 
of R and p had no effect on the rate of convergence 
of the eigenvalue spectrum. For the angular 
expansions in confluent hypergeom. functions the rate 
of convergence is by far inferior to the one found 
using expansions in ultraspherical polynomials. The 
table listing the number of significant figures of the 
separation constant generated by a 15-term expansion in 


confluent hypergeom. functions 


Sw(x) % q 2 (c3,n;x) [226] 
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The cases where the number of significant figures is 
indicated as ">16" are the ones for which the expansion 
¥sttruncated) (at. less than 15 terms) due to a coefficient 
in the recurrence relation being zero. 

If the first parameter of the confluent hypergeometric 
functions serves as expansion index the radial expansions 
are obtained. There are eight of the 16 expansions proposed 
Pie 2 6.) sot thisetyper Themsitiation 31 sais) larto 
the one observed for the angular expansions in confluent 
hypergeometric functions. There is a total of four distinct 
eigenvalue spectra of the truncated radial expansions. Only 
two of them converge upon extending the number of terms 
considered in the expansion. All four of them can be 
generated taking one block of an arbitrarily chosen 
expansion in confluent hypergeometric functions considering 
themcour DOsSslore COMOINatIOnSs 1h. slgn LOL = Ro andy. Dp. 

Unlike the angular expansions the finally converged 
eigenvalue spectra (the few lowest eigenvalues !) are 
distinct. The eigenvalue spectra of interest are the ones 
where the wave function remains bounded for \>*+@. 

This is, the case if for the subblock for the polynomial 
solutions of the confluent hypergeometric DE the signs 


are chosen as follows 
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tie ton a gaven-combination of signs of R and 7, 
the sign of p is changed then the eigenvalue spectrum 
of the infinite sub-block located opposite to the one 
considered in the first case will be identical to the 
original one. Hence if non-polynomial solutions of the 
confluent hypergeom. DE are employed in the expansion 
i.e. n2l the proper combination of signs will be the 


one obtained by changing the sign of p in [227] 
eel, © oe sh Be [228] 


The rate of convergence for this expansion is the 
same as the one of the Hylleraas and Jaffé expansion, 
although the generating function used by Hylleraas 
is not exactly the one of the Laguerre polynomials. The 
equivalence of his own expansion with a Laguerre 
polynomial expansion was apparently already known to 
Hylleraas. Four tables are given below listing the 
number of significant figures of the separation constant 
for different combinations of parameters of the GSW-DE. 
As was expected on the basis of qualitative speculations 
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Tables of the actual numerical values of separation 


constants for the cases listed are provided in appendix II. 


Since Laguerre polynomials are related to Hermite 
polynomials as follows 
ale 
Hy (x) = (-1)°P (ntl) Lo 4 (x?) [229] 


Wea eet 
Done eae y 


L 
T (nt+1) xLy (x?) [230] 
one might suspect that expansions in terms of Hermite 
polynomials and functions would generate the radial 


solutions to the GSW-DE. 


Looking for expansions 


y gq) ¢ (x2-p2) D242 (m+1) xD-x2+m (m+1)-A} Ho i) 60, A237) 
rire att 3 
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£23.2)] 
using the relations, 
x°D*H (x) = n(n-1) (H, (x)+2 (2n-3)H,_. (x) 
+4 (n-3) 5H, (x) ] [233] 
2 
xD HL) = EE pelea eel ae Seren) as e0Ge ls [234] 
eis) = aera eS [235] 
2 y ul 
x DH (x) = n{(5)H 4) (x) + (2n-1) 8, _, (x) 
+2 (n-2) 5H _ 3 (x) ] [236] 
xDH_ (x) = nH. (x) +2 (n-1)H,_. (x) J e234 
DH, (x) = 2nH_, (x) [238] 
poncho gee 
ed aie wt (0) tint ()) (oc) 
+(n-1) 5H _» (x) [239] 
ih 
xH (x) = (S) HL 4 (x) tnH_y Cs) [240] 


one readily obtains recurrence relations for the 
expansion coefficients in Poo meander 2 oe. 


(1) 


n 


CL) 


(1) + [m(2n¢m+1) +n? ~ (5) -Alg 


ul 
rene 
+[(nt1) 5 (4 (ntm-p?) +5) 1g P+ 14 (nt) plant = 0 [241] 


(7) £05) (kt) J4q,07) [n (ntby-1) kg] tay2) [ (n+) 


(2 (bytn(a,-1)) +k, -1) 144,62) [2 (nt1) 9 (2n+b, +1) ] 


+q(?) [2 (n+l) ,(2a,-1) I+a\f [4 (ntl) 4] = 0 [242] 
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tn [242] thessubstitutionspli sted) inwl221]=[217)] 


would have to be made. The eigenvalue spectra derived 
from [241] and [242] are the ones of the angular 
eignefunctions. For moderate p* they were found to be 
very slowly convergent. Using a 15-term expansion in [231] 
hor-diffterent QO,'s and p's it was found in general that 
more than half the eigenvalues of the obtained spectra 
were complex conjugate pairs. This constitutes further 
evidence for slow convergence of these expansions 

making them entirely unattractive for the use in 


numerical work. 


From Fig. II it is apparent that the expansion in 
[226] converges faster with increasing internuclear 
distance which indicates a possible usefulness of the 
angular expansions in confluent hypergeom. functions 
if corresponding sections of the potential energy 


curve of a heteronuclear species are to be explored. 
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Jaffé Expansions for the GSW-DE 


A bilinear transform like the one used by Jaffé can 


be applied to any GSW-DE : 


{ (x?+a,x)D*+(bo+b xt+b,x* ) D+k +k, x} y(x) = 0 [243] 


ih 0 
This transform is achieved in three individual steps 
as shown below. First we divide out an arbitrary power 


‘ome, pe ae 
TAK) ke Vi) [244] 


which transforms [243] into : 


{x? (x+a,)D°+x[2aa +by+ (b, +20) x+b er Dee tob,)x* 


ji 2 


i 
+ [a (b, +a-1)+k 9] x+al (a-1)a,+b)] } (=O 


[245] 
Then a is chosen such that the quadratic term in 
the coefficient of j in [245] vanishes 


a = -k, /b, [246] 
Substituting for the variable 
x =t [247] 
[245] transforms into 


(ee (tat Desh 2ieol ja tb) te (b te 2Casl)esbo |p 


+[a((a-1)a,+b Jt+a(b,+a-1) +k} inte) = 0 [248] 


0) 
Making the additional substitution 


t= (s-1)/a, [249] 
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the last DE transforms into 


Ve (eal) oD] (2 (ol) a sembyn2(Tto-1) )sta bo beer 1D 
ta[T)+a-1] sta(b)-T5)+ky} j(s) = 0 


where To = b/ a, [250] 


A power series expansion for j 


j(s) = 2 ae P2515 


is defined by the three-term recurrence relation 


among its expansion coefficients 


Gy b (n-1) (n-2a-T)) +a (Tota-1) ]+q_ [n(2 (T5+a-n) -b, ) 


ta(b)-T))+k I +q [(nt+1) (ntb, -a,b5-To) ] = 0 WA 


Tote 
which was obtained by substituting [251] into [250] and 
applying standard procedures. Notice that the indicial 


equation of [250] has roots 


=O, CaS T,+a b.-b +1 [253] 


=i Omer i 


If for [250] we assume a Frobenius type solution 
eS) Sena ees = s’h(s) [254] 
[250] transforms into 
Pees) eDet sl 2 (batty ee 2(T hae boi) Sba) S 
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For h we can have power series expandable solutions 


ns 
By = yi) 
ae B. = Tota, b,-b) +1 
B, = a+T, 
By = atl [256] 


and they are defined by the following recurrence 


relation among their expansion coefficients : 


dy-] [n(28-a, b,+b, tn-1-T)) +8 (B-a,b,+b, -T)-1) ] 


tq [n(2(T +ta-2B-n) -b, ) +a (by -T))-8 (2(8-T,-a) +b, ) +ko] 


0 
a2 Gras [n(n+1-T)+2 (B-a) )+8 (B-20+1-T)) +a (LT yta-1) ] = 0 
E2574 


Feasibility Study for the Frobenius-Jaffé 


Type Expansions 


Examining the limiting form of the radial expansions 
for R+0 presented in [29] we realize that the radial 
wave functions might asymptotically approach non- 
integer powers in the neighbourhood of the expansion 
center for small R. 

In this situation rapid rate of convergence ican be 
expected to occur if at the expansion center the 
expansion terms (or at least one of them) show the same 


asymptotic behaviour. 
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Jaffe did not mention that his DE had an indicial 
equation with roots differing by a non-integer number. 
Apparently he rejected Frobenius type solutions for 
physical reasons, i.e., the shape of the radial wave 
FuNCe ONnseatarnesnucter i orem—0. 

It can be shown that a non-integer power at the 
expansion center can be generated with a Jaffé type 
expansion. In this case it would do no harm in respect 
to physical acceptability if a non-integer power is 
factored out. To show what is meant by the above statement 


we can examine the second order equidimensional DE 
{x? (a? /dx*)+ax (d/dx)+b} g(x) = 0 [258] 
havingvarchar. polynomial with the pair vor Toots 
ele (ena) ty (ae) eb [259] 
Substituting for the variable 
s = x/(x+2) [260] 
transforms [258] into 


(e* (s=1)- Hs (S-1)" 2e-ayDrbr eo (s)) = 0 [261] 


a DE which has an indicial equation with the same 
LOOts asi thewindicual equation of [258 7.2 4e., the 
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A power series expansion for g 


g(s) = z qn oe [262] 


is defined by the following recurrence relation : 
Pata = 2 58, 
q,-1 [7 1] gq, [n(2(n+1)+a)]+q_,, [n +n(2+a)+(2+b)] = 0 
l263)] 
Notice thatein the tradiag. (tridiagonal), matrix 
corresponding to this recurrence relation the diag. 


20 


elements grow with "-2n whereas the upper and lower 


off—-diag. elements grow with “n?". The same is the 
case for the radial expansion in confluent hypergeom. 
functions and the Jaffé type expansions. 

This is an indication that these expansions actually 
do build up a non-integer power at the expansion 
Cencerr. 


Notice as well that non-integer powers and logarithmic 


functions can be expanded in confluent hypergeometric 


functions 
a : a 
x = DTtateth) F 1 (-a)_/T latnt)) Lot) 
n n n 
O<n<~m , -a<l+min(a,a-%) [264] 


Trey = eee atl) zr, [T (n) /T (atnt1) }L2 (x) [265] 
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Jaffé Expansions of the Original GSW-DE 


If in [256] we make the substitutions [212]-[216] and 


[250] and taking bj=-1 we get : 


By = 0 [266] 
B. = -m [267] 
B3 = (B/p)+m+1 [268] 
By = (B/p)+l [269] 


The pairs B85 and B4/8, dit fter sby ee Maawi ches 
an integer. Hence the expansion of [255] with B=B, is 
equivalent to the one with B=B., if we just replace 
the expansion index "n" by "n-m" and the same relationship 
applies to the pair B3,8,- 

In the original GSW-DE as for example written in [18] 
replacing "R" by "-R" is equivalent to changing the sign 
of the variable. From the form of the Jaffé expansion as 
shown in [48] it can easily be seen that taking B=atl 
(or B,=o+1 if referred to [48]) is equivalent to changing 
the sign of "R" and taking expansions with B=B, or 
equivalently B=B.. 

Hence due to the six coefficients in [243] being 
interrelated by only three parameters R,m and E, we have 
actually no additional option if we resort to Frobenius 


expansions of [255]. 
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The table below summarizes this situation. The 


eigenvalue spectra are designated with {A}, with 


Poorer 4. 

k Peel eee. iy 
p=-|p| , R=-|R| 1 1 gt at 
p=-|p]| , R= |[R| 2 2 1° rig 
p= |p| , R=-[R| 3 3 4t me 
p= |p| , R= [RI 4 4 go gS 


{A} 


1 is the eigenvalue spectrum of the physically 


meaningful eigenfunctions of the radial GSW-DE. 

The matrices arising from B3 and By are the 
transpose of the ones obtained using By and B.- This 
is indicated by the "f‘s". This notation holds 


throughout. 
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CHAPTER ITI 


Integral Transforms 


It is a well-known fact that integral transforms 
have played an important role in the investigations 
and solutions of many different types of DE's. 

In the present case the wave function would simply 
be the preimage of the expansion of the transformed 
GSW-DE. One would expect that if we guarantee the 
existence of the integral transform by formulating 
an algebraic eigenvalue problem as before, in the 
case of convergence for a given set of parameters 
R,p and m, we guarantee simultaneously the existence 
of the preimage of this expansion. In other words, 
the same eigenvalue spectra are expected to occur 
when working with integral transforms as by 


expanding the original wave function. 
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The Laplace Transform of the GSW-DE 


It can be shown that by a Laplace transform 


a GSW-DE is transformed to a GSW=-DE. 


We are interested in the following transformation 


£{S(L(A)]|s} = £{[ (A?-1)D?+2A (m+1)D- (pa)? 


+20QA+m(m+1)-A]L(A)|s} = 0 [270] 


where according to the distinction between an 
angular and a radial expansion the following 


ambiguity of sign occurs : 
Qi 4(2Z,4Z5)R [av 


The Laplace transform of the wave function 
which we write as £(s) satisfies the Laplace 
transformed DE in [270]. The transforms of the 


individual terms are 


£{(X2-1)D*L(\) |s} 
= (d?/ds”) [s*2(s)-sL(0)-L' (0) ]-s*2%(s)+sL(0)+L' (0) 
[oe 


£{2A (m+1)DL(A)|s} = -2(m+1) (d/ds) [s2£(s)-L(0)] 
ena 
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£{[-(pdr) *+2Q\+m(m+1) -A] L(A) | s} 


= [-p? (d*/ds?)-20(d/ds)+(m(m+l)-A)]2(s) [274] 
The transformed GSW-DE then becomes 


{ (s*-p* )D?-2 ( (m-1) s+Q) D+ [-s?+m(m-1) -A] }2 (s) 
= -[sL(0)+L' (0) ] [275] 
for which we seek the homogeneous solution. 
Applying the transformations 
ees) a exp(s )g, (s) te2 6) 
Ls) = exp (-s)g, (s) [207 


one obtains a pair of DE's 


{ (s*-p?)D?-2[s?+ (m-1) s+ (Q-p*)]D 


+[2(m-1) s+m(m-1)+20-p?-a] } sel = © [278] 


{ (s*-p?)D*+2[s*-(m-1)s-(Qtp*)]D 


+[2(1-m) stm(m-1)-20-p?-A] } gS) = 0 [279] 


Shifting the center of expansion 
ase = Bae [280] 
where we choose the plus sign. This available choice 
is equivalent to choosing the sign of p, an option 


available to us for the original GSW-DE as well. 


{x (x-4p)D?-[x?7+2 (m-1-2p)x+4 (Q- (m-1)p) JD 


+{ (m-1)x+(m-1) (m-2p)+2Q-p*-A] } cies) = w [281] 


{x (x-4p)D?+[x?-2 (m-1+2p)x-4 (Q-(m-1)p) JD 


+[ (1-m)x+(m-1) (m+2p)-20-p*-A] } g(x) = 0 [282] 
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The last two DE's are equivalent : [282] can be 
obtained from [281] by changing the signs of p and R and 
the variable x. A power series expansion for Ty 


gy (x) = Boqux" (2834 


is defined by the three-term recurrence relation 


holding among the qn s 


qj] m-n] +q, [n (n-2m+4p) + (m-2p) (m-1)+2Q-p’-A] 


= 4q_,,[ (ntl) (Qtp (n-mt1) J [284] 


The expansion terms can be back-transformed if for 
the expansion index the following is specified : 
-l<n<-« [285] 
which gives 


SOS pl see pies [286] 


For inverse Laplace transforms in general the 


LOLIOWING yholds =: 


£*{a,f, (s)ta5f,(s) |x} =a £{£,(s) |x}+a,¢ 7 {£,(s) |x} 


i 
(23,7) 
However if the sum has infinitely many terms this 
procedure, being a combination of two limiting 
processes, becomes problematic. Nevertheless it 
is believed that, if the eigenvalue spectra of the 


expansion of the Laplace transformed wave functions 


converge, then existence of the inverse Laplace 
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vist 
transforms of all expansion terms is sufficient to 
assure that what is obtained by term-by-term back- 
transformation represents the Laplace preimage of this 
expansion. 

A nice example where we find a rapidly converging 
angular eigenvalue spectrum but cannot achieve 
term—by-term back-transformation is an expansion in 
Jacobi polynomials shown below. 

FOR dnwexpanstonmeote [27 Sise(ands (27.90) )esins Jeacobr 


polynomials the following relations are needed 


(2nta+B) xP 98) (x) = 2 (nt1) (ntatB+1) (2ntate)P (4%) (x) 
-(2ntatB+1) (a?-B2)P 1") (x) -2 (nta) (nt 8) 
(Qnto+6+2)P (418) (x) [288] 
(2n+atB) (1-x?) (ayaxyen ce (x) = [a-B-(2ntat+B)x]P\%" B) (x) 
+2 (n+) (ntg)P 178) (x) [289] 
and the linear differential operator satisfied by 
the Jacobi polynomials , 
{ (L-x?)D?4 [B-a- (atB+2)x]Dtn (ntatp+1)} PL%"®) (x) = 0 
[290] 
The parameters a and § were chosen such that 
bothersome terms vanished as follows : 
oe = "=(m+(0/p))) 2945] 
B = -(m-(Q/p)) [292] 
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The expansion , 


g,(s/p) = 2 a.Pi%’®) (s/p) [293] 


if substituted into [278] generates the following 


recurrence relation for the om Siar 


,-1 [2np (2m-n) / (2 (n-m) -1) ] +g, [n (n-2m+1) +m (m-1) -p?-A] 


See 2b (O/ py moma t) 7 (2 nem) 3) [294] 


The truncated characteristic polynomials derived 
from [294] are the same as for the expansion in 
ultraspherical polynomials which themselves are special 


cases of Jacobi polynomials, 


(at (4)) 
Cr 2°’ (x) = [[T (atl)? (2atn+1)]/[T (2a+1)T (atn+1) ]] 


pieert(x) [295] 


The Jacobi polynomials do not have inverse Laplace 
transforms. If one chooses, 

-l<n<-« [296] 

the Jacobi functions can be formulated as generalized 


hypergeometric functions : 


ay (x) = [[ (atl) /P (ntl) ] FP) (-n,otptntlsotl; 


(l-x) /2) [297] 


If we make use of the well-known linear transformation 


formula, 
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and using the formula 
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recat yr (oy) 


fag Iam SHea(ary amribane/ac)ilitd 

oF 5 (ay 1 agib, . 07 2t) : Re (oa) >0 [29 | 

and the shift theorem of the Laplace transform,one 
encounters further difficulties since,for the second 


fermeinal2clesthe*restoietion for oUin | (299 lmdoes not 


nNoldsisinvisechosenvas /shownkin (296). 


ne roots 00 thie indicia l equationssor @|/238) mand 
[282] in general differ by a non-integer number. 
They are 
Coes Qe te MS (0/D} [300] 


We substitute in [281] , 
Hes) = ot On es [301] 


where the substitution "“T=0/p" is made throughout 
this section. 


Then [281] becomes , 


{x (x-4p)D?-[x*+2 (T-2p-1) x-4 (Q-p (m+1) ) ]D+[ (T-1) 


(T+x-2p)+2mp-p*-A]} h(x) = 0 [302] 
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A power series expansion of h, 
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is defined by the following three-term recurrence 


relation : 


Apps llabsied estes Veal asks 9s) a vinlely | auesy yey l(a eel.) 


+2q,4 , 0 (ntl) (Q-p (ntmt+1) ) ] = 0 [304] 


Two distinct eigenvalue spectra can be obtained 
from the last expression depending on the choice of 
the signs of R and p. The same eigenvalue spectra are 


obtained if from the original GSW-DE, 


{x(x-2)D*-2 [px* - (2ptm+1) x+ (m+1) ]D+2[Q-p (m+1) ]x 


+[m(m+1)-2(Q-p(m+1) ) -p?-A] } Geax lia 0 [305] 


a power is factored out , 


x (T-m-1) 


93(x) = hy (x) [306] 


which gives , 


{x? (x-2)D?-2 [px?- (T+2p) x+ (2T-m-1) ] xD+ [ (T-m-1) 


(T+m+2p)+m(m+1)-p*-A]x-2(T-1) (T-m-1) } h, (x) = 0 
[307] 


A power series expansion of hy " 


i n 
h, (x) = es qn* [308] 


is defined by the following recurrence relation 


-2[(n-1)p]q,_,+{n[2(T+2p)+n-1] + (T-m-1) (T+m+2p) 


+m(m+1)-p?-A]q-2[ (n+l) (2T+n-m-1)+ (T-1) (T-m-1)]q, 44 
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One eigenvalue spectrum of [309] is obtained if the 
signs for p and Q are the same and a distinct one is 
found if these signs differ. The two eigenvalue spectra 
although being convergent apparently are physically 
meaningless. 

Since a Laplace transform generates a GSW-DE from a 
GSW-DE,one is tempted to explore “higher generations" 


of integral transforms. For example writing [302] , 


(Gla Cleve Obie )ahen PeNGu tio) oal i (esvies sh) Woe" te 


+[ (T-1)/x]+[[ (T-1) (T-2p)+2mp-p*-A] /x*]} h(x) = 0 


[310] 
and taking an inverse Laplace transform , 


£71 {n(x)|s} = 2(s) 
gives the following integral equation : 


{s(stl)2(s) - ih [4pt2-2(T-1-2p)t-(T-1)]2(t) at 
- S65 fg (4(Q-p (mt) ) t+ (T-1) (T-2p)+2mp-p?-A] 


‘bite trek pipe aa fenee 
Taking twice the derivative gives the DE : 


{s (st+1)D?*+[T+1+2s (T+1-2p) -4ps”]D-4 (Q-p(m-1))s 
+(T+1) (T-2p)+2mp-p*-A} 2(s) = 0 [3i2] 
If from the original GSW-DE one takes twice the 


Laplace transform the original GSW-DE is obtained 
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with the sign of the variable changed. 


If in [270] we make the substitutions 


t = exp(ah) (313) 
(d/da) =sat(da/at [314] 
(d/dn) * 9a" (fo (ad /dt) 7+ (dZar)) [315] 


we are dealing with equidimensional linear 
differential operators. Then making use of the 


following two general properties of the Mellin 


transform : 


{ (x (d/dx))o (x) |z} = (-z)™M{o (x) | z} 


[316] 
M{(1n(x)) "9 (x) |z} = (d/dz) "M{¢ (x) |z} 
(3187, 
we can generate a second order linear DE with 
polynomials as coefficients from [270]. 
The Mellin transform of ¢ iS given by 
co 2-1 
M{o(x)|z} = fox (x) dx [318] 
and it is related to the Laplace transform, 
which is apparent upon performing the simple 
pete : 
substitution aL ae, Cages [319] 
which gives a two-sided Laplace transform 
M{g(x)|z} = S.2(e “ye at [320] 


or equivalently a sum of two one-sided Laplace 
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integrals of parameter z and -z 


co 


0 


EVA 


Mid: Cail izigeeerye Me oe dt 


i pi eNeer ae [321] 


After performing thecusubstitutions Listed in 


[313]-{315] the limiting factor of the Mellin preimage 


will be 
T(t) 9 teRAe) | t | + [2223 


and [270] becomes 


{[1n* (t)-a?] (tD) ?7+2 1n (t) (m+1) (tD) 
~ (p/a) *1n* (t) +2 (Q/a) 1n(t)+m(m+1)-A} L(t) = 0 
i323) 
Denoting the Mellin transform of the wave 
function 


M{L(t)|s} = g(s) [324] 
the transform satisfies 


{ (s*=-p*)D*-2[ (m-1) s-Q]D-s*+m(m+1)-A} g(s/a) = 0 
[325] 


which is the same as the Laplace transform 


[275] except for the sign of Q. 
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Expansion of the Laplace Image in Confluent 


Hypergeometric Functions 


In the previous section it was shown that the 
Laplace transformed GSW-DE is itself a GSW-DE. Since 
the expansion in confluent hypergeom. functions were 
worked out for a GSW-DE with unspecified coefficients, 
the recurrence relations [196]-[211] can be applied. 

First the inverse Laplace transforms of the 
expansion terms will be examined. We need essentially 
four back-transforms to cover all the expansions 
proposed, if the original Laplace transformed GSW-DE's 
[281] and [282] are expanded. 

Expanding [302], where a power corresponding to the 
non-zero root of the indicial equation was factored out, 
two additional inverse Laplace transforms have to be 
evaluated. The exponential factors occuring in the 
expansions are taken care of by simply employing the 
shift theorem of the inverse Laplace transform. 

From the literature by standard mathematical operations 


Sie we extract readily the following result : 


Pane (gyre) tale oe tele ey (oe Cee 


Re(a)>O , Re(s)>0 [326] 
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If the exponential factor of the expansion is taken 
into account we find that term-by-term back- 
transformation of the expansion according to [326] the 
familiar Jaffé expansion is obtained. From this we gain 
some understanding regarding the initially perplexing 
equivalence of the convergence rates of these entirely 
different expansions. 

The remaining three inverse Laplace transforms are : 


1 


£"(8(o,yis)|t} = t'5F, (a,lryst +) [327] 


atts y(a,yzs) fe} = te sr (a-yth) 1, Fo (ar et) 


af (a-y) A 2 ee 
Sed re (Gs ytL)] Zyl (a) 71) I i) 


Poteet Gea WinALL SeeRelioey 41) SO eRe seo 
[328] 
st ts) 6 (a, y;s) |t} = [P(y)/t"] [t-1] % [329] 


The last two formulas indicate again a Jaffé 
expansion as the Laplace preimage of an expansion in 
confluent hypergeom. functions. 

The inverse Laplace transforms above were derived 
from tables of Laplace transforms by Oberhettinger 


and eae 


making proper substitutions. 
Finally we consider the two additional cases 


arising from confluent hypergeom. expansions of [302]. 
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The inverse Laplace transforms [226]-[229] are 
special cases of [330] and [334]. In the tables 
listed below, the eigenvalue spectra of the 
eigenfunctions of [281] for the different expansions 
and sign combinations are shown. 

For the radial expansions, the numbering adopted 
is the same as was used in the Jaffé expansion of the 
original GSW-DE. 

The parameter "k" refers to the expansion Sy used 
and the numbers in the tables below are identification 
tags for the eigenvalue spectra {A}, , where {A}, is 


again the physically meaningful radial eigenvalue 


spectrum. 
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The numbers on the left of the commas apply to the 


infinite subblocks, where for the expansion index n 
we have n>-~ and the number on the right holds for 
the infinite subblocks corresponding to n>+-. 

If m#0 eight differently converging angular 


eigenvalue spectra are found. They interrelate as 


follows 
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For the Y-states the following special situation 
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In the table above the sign combinations were coded 


asmshowectn | [4//3255 | as 


wien ths 


ul 
D 


Jaffé Expansion of the Laplace Image 


Since the inverse Laplace transform of the Jaffe 


expansion is 


Per teeet (s-1) Ast) eh = ae 


exe O(a) | (335) 


one would suspect that the same rate of convergence 
is obtained as for the expansions in confluent 
hypergeometric functions. It was indeed found that 
the previously found four radial eigenvalue spectra 
are generated. It appears that integral transform 
techniques provide no cure for the convergence 
problems encountered for the radial expansions. 

For the different Laplace transformed GSW-DE's 
the algebraic eigenvalue problems associated with 
the Frobenius expansions in [256] interrelate as 


follows : 


Expansion of [281] : 
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Expansion of [282] : 


iAyL By Bo B3 By 
p=|p| , R=|RI 4 3¢ 3! 4 
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Expansion of [302] =: 
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The {A}, 's in the above three tables correspond 
to the {A}, 's in the table following [269]. 

In the cases where a superscript "°" is 
indicated for the 8's the expansion index n has to 
be replaced by "n-m" or its range has to be changed 
accordingly. 

If m=0 the matrices for the Laplace image expansions 
are real symmetric. This is the case as well for [312] 
where an additional Laplace transform was taken. For 
this DE the truncated radial eigenvalue spectra are 


interrelated as follows : 
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Expansion of [312] : 


{A}, BS Bo B3 By 
p=|p| , R=[R| 3 Aa 37 4 
p=|p| , R=-|R| 4 37 at 3 
p=-|p| , R=|R| 1 gt 1" 2 
p=-|p| , R=-|R| 2 2t at 1 


However for m+#0 the off-diagonal elements for 
these are different from the ones of the Jaffé 
expansions of [281],[282] and [302]. In the cases 
where a "t" is indicated the corresponding matrices 
are adjoint to the ones where this "t" is missing. 

The Laplace transform acts as the Rosetta stone 
to interconnect the expansions in confluent 
hypergeometric functions and Jaffé's Mobius 
transform. However no change in the convergence 
properties of the separation constant results 


from the use of Laplace transforms. 
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CHAPTER IV 


In this chapter results of calculations on the 
particle in an electric dipole field of two point 
charges are given. 

Rather frequently people consider free electrons 
in the company of polar molecules as electrons 
interacting with dipoles. Two examples of 
idealizations of physical situations resulting in 
the treatment of this system are : electron capture 
by rotational excitation of polar molecules ew) 
theory of low-energy-electron scattering by polar 
molecules, ‘1+) 

It was known previously that the binding energy 
Of a particle in a dipole field becomes zero at 


oy Lithia 


some critical internuclear distance. 
chapter the "dipole theorem" is presented. It 

allows one to readily evaluate the critical distance 
for any combination of quantum numbers. Although no 


proof of this theorem acceptable in every respect 


will be shown, the correctness of this theorem is 
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considered to be established beyond any reasonable 
doubt. No extensive set of exact energies for this 
system existed previously. Results for ten states 
of the particle in a dipole field are presented 
here. As thes critical lpointsr on tie porencial 
curves are approached numerical break-down of the 
radial expansion is encountered. The program 
performs a test for final convergence by extending 
the number of expansion terms considered and 
comparing with the previous result. Although the 
algorithms used in the computer program allow one 
to include essentially any number of terms in the 
expansions, it is believed that 45 expansion terms 
constitute an upper limit beyond which convergence 


tests such as the one mentioned become questionable. 
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The Particle in a Dipole Field 


Le 12, |=12,| , we are dealing with a special case 
in the sense that for a given set of quantum numbers, 
the potential curve is related to the potential curve 
of the same energy state belonging to the pair of 


charges |2,|=|2 with 


4| 
[336] 


as shown by the following operations : 


To obtain the energy belonging to Zs and Zo 


at Ry 2 take the energy belonging to Z. and 
g 
VAGE Gebel 186 where 


R34 = Ry 2/t [a3 7.) 


and divide this energy by the square of the 
Sea iting Fractor ss 


= 2 
E = E3 4/5 [338] 


Lee 
ile Z,=-2, we have the situation discussed above. 


The radial GSW-DE then becomes 
{ (A?-1) D*+2 (m+1) AD- (pA) *+m(m+1)-A, } Los" 0 [339] 


Inspection of the above DE readily indicates that 
as p>0 the L's exhibit the following asymptotic 


behaviour : 
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where the Z's are modified Bessel functions if 
the limit is approached from the side E,” 0- 


TEethe substi tution, 
Se [341] 
Vs made in [339] the following DE results 


{t (t=1) D241 (m+ (3/2) )t-(1/2) 1D=(p/2).74 
+{(m(m+1)-A,)/4]} L(t) = 0 [342] 


Since one root of the indicial equation of 
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[342] is a corresponding power can be factored 


out 
lines) sg tectar te) [343] 


which gives 
PeeWee 5/2) (3/2) Oo py 2) ot ew 1/2) 
+[ (m(m+3)-A,)/4]} h(t) = 0 [344] 
Power series expansions of L and h such as 


ie) Sens eee [345] 
n=0 = 


ear) 4a Ce ac Oe [346] 
n= U4 


II 


TaG) 
give rise to the following pair of recurrence 


relations a 


(1) 


-p?q 1) +[4n (ntm+ (1/2) )+m(m+1) -A,] a, 


-4[ (n+ (1/2)) (ntl) 144) = 0 [347] 
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(2) 


n-l 


(2) 


-p*g +[4n (ntm+ (3/2) )+2+m(m+3)-A,]q, 


-4{ (n#1) (nt (3/2))]q7) = 0 [348] 


As expected the eigenvalue spectra arising from 
[347] and [348] are the angular eigenvalue spectra 
for the bonding and antibonding states. 

A rather interesting aspect of the use of Laplace 
transforms for the spheroidal wave functions is given 
by the following 


If one takes the Laplace transforms, 


Site) Ste. S) [349] 


Sitti) | Sie ys) [350] 
where the transforms satisfy, 
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then the power series expansions 
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are defined by the three-term recurrence relations 
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whose eigenvalue spectra are the ones for 
antibonding states if the Laplace preimages were 
bonding wave functions and vice versa. 

This phenomenon can be explained by term-by-term 
inversion of the Laplace image, 
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and the relation of the Laplace transformed 
GSW-DE to its preimage for the case Q=0. 
If p=0, [342] (also [344]!) becomes a hypergeometric 


DE2BwOnemsolutioneis 


L(t) = (m+ (1/2)+ (A, -(17/4)) 977) 72, 


gig (A 

((m+ (1/2) -(A,-(74)) 27?)) 72; 

sal 2) te) [358] 
which is however unconditionally divergent if m20 

Onethesboundany Of (thevunitg disk yan “the complex Pplane. 
The limiting form of the radial wave functions jas 

indicated in [340] suggests the use of a Neumann series 


expansion. 
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Fa (x) [361] 
[360] transforms into 


{ (x*—=p*)D*+ [xt ( (2m+1) /x) p? ]D-x?- (1/4) -A, 


-(p/ (2x) )* (2m+1) (2m+3)} Fa(x) = 0 [362] 


then using the following Neumann expansions : 
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one obtains the following three-term recurrence relations 


among the expansion coefficients in [363]-[366] 
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where the minus sign refers to the expansions in 
modified Bessel functions and the plus sign refers to 
expansions in Hankel or common Bessel functions. 
If modified spherical Bessel functions are 
employed with 
Jn | 2m+5 for bonding states 
|n| 2m+3 for antibonding states 
then the eigenvalue spectra for the angular expansions 
are obtained. The rate of convergence of these 
expansions is rather impressive. Expansions of the 
second lowest angular Y-state with four and five 
expansion terms gave for the separation constant 


the number of significant figures shown in;the 


following table 
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Dp 0 5 10 15 20 25 30 35 40 
4 terms >16 9 8 7 6 6 5 5 5 
5 terms >16 13 1! 10 9 8 u 7 i 


Bromeinvestigations: which will be shownslatersin this 
section, evidence was obtained indicating that the 
radial expansions would be obtained if use could be 
mMmeadesor Bessel srunctions Of integer order. This can 
be done trivially if p=0. Furthermore such an expansion 
is possible if m=-5, The latter case, although 
physically meaningless is in a mathematical sense of 


distinguished importance, since for this case a 


subst butlon 
AX = cos (z) [376] 


in [339] gives a Hill-type DE known as Mathieu's 
Dia Lt me-> the expansion in Bessel functions of 
integer order cannot be obtained due to occurrences 
of inverse powers of the expansion index in the 
recurrence formula. 

Later it will be shown that as p-0 the radial 
eigenvalue spectrum for the separation constant 
collapses into one infinitely degenerate eigenvalue, 


thus explaining the farlure of the present approach 
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because the DE satisfied by Z (Pr) becomes an 
equidimensional or Eulerian DE whose characteristic 
polynomial has the roots tv. 

Therefore, if for p=0 we assume a trial function 


a Ae 
y Nels) eee [378] 


then for a proper exponent e€ the correct 
eigenvalue spectrum should he obtained. 
Taking the Rayleigh coefficient as a first 
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AP 
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/$°x"?£ an] = ef ((etl) /(e+5))- 
See e>t [379] 
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Lr O.N)), 7a (q,tq,A )A [PS e2e5) 
the following Ritz-Galerkin eigenvalue problem 
arises : 
I i eri =(e+2k) 
iZo non C5;Cy, <h [sa Sip Ones S| 


with according to the Hylleraas-Undheim theorem 


improved eigenvalues. As before the differential 


operator of the spheroidal wave functions is written 
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Square: integrable, this valine is still) acceptable nO. 
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Subst Peet onsodk c=5 in the expression for the 


Rayleigh coefficient gives 
A, = m*-— R39 


Another less laborious approach for finding the 
characteristic exponent e« is based on a theorem of 
upper and lower bounds of the eigenvalue established 
too 2B by ery orm sand Bogoliubov '?). We state this 


theorem as follows : 
Having an eigenvalue problem 
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According to the theorem of Kryloff and Bogoliubov 
the interval in which at least one eigenvalue is 
located has zero width if a*=@*. In order to find the 
values for e€ where this holds one would have to locate 
the roots of some rational function. However based on 
the previous investigation we are particularly 
interested in finding out what happens if e=5. For this 
case a=+f becomes 
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where obviously the plus sign holds. 


These investigations seem sufficient to establish 


the following theorem 
Theorem : "Dipole Theorem" 


Vie E ,=0 the separation constant 
eigenvalue spectrum of the radial 
expansion of the spheroidal wave DE 


becomes infinite-fold degenerate with 
ree ae [398] 
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radial quantum number. This result, although ina 
mathematical sense arising trivially, is rather 
unexpected in a physical sense because one would 
expect the binding energy of the particle to become 
zero if the Bohr radius of the hydrogenic system is 
appreciably penetrated by the additional charge. 

The internuclear distance where the binding energy 
OGraneeclectronrin a o1pole t1eldeis  .7eroutseot 
considerable importance to radiation chemists working 
with negatively charged “heavy" particles. For a 
negatively charged "heavy" particle in dense matter, 
this distance gives us an indication with respect to 
the diameter of the tunnel along the trajectory of the 
Damticle@in whose. interior JOnization of vardouss.types:. 
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OLD electrons can occur. Turner and Fox 
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function used by these authors is very interesting. 

The author calculated the value OF at p=0 for ten 
angular states where use was made of the previously 
stated "dipole theorem" and the angular expansion in 
uUbEeraspherdjcal polynomlais. Ihese data, are presented on 


the next page. Turner and Fox indicate their result with 


Six figures. However the last “two of these are incorrect. 
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|Z, | ie |Z, | = (399) 


1.e.,a proton and a negatively charged meson, 
ae is the actual internuclear distance, Otherwise 


the formula 


ea 2 eee RS OLIVA) On Zn [400] 


has to be used. 
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qualitatively would be expected since for large m's 
the magnitude of the wave function in the region 
along the internuclear axis and hence in the 
neighbourhood of the negatively charged meson is 


small. 


An extensive set of data was obtained by the 
author for the states 32n>2%2m. The first point in 
the corresponding ten tables was obtained by simply 
applying the "dipole theorem". We recall that this 
theorem was not established on a firm basis because 
the given eigenfunction was not square integrable 
over the physically meaningful domain. However the 
points derived from the "dipole theorem” fit nicely 


into the calculated energy curves. 
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One would think that the energy of the system 


rapidly rises as soon as the negatively charged meson 


PeNnetrates tnesbOniraraciiseoO1 thee leccron .siOwevier 


the data above indicates that the energy rises gradually 


in such a process. 


The ten dipole energy states which have been 


treated were identified as follows 


oo 
First a point on the E,(R) curve for the H, |n; &,m > 


energy state close to the united atom limit was 
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evaluated. Then as points for larger internuclear 
distances were being calculated a small increment was 


simultaneously subtracted from Z, until Z,=-1. This 


2 2 

gave one single point for a large internuclear distance 

on the EB, (8) curve of the |n;2,m> state of the dipole. 
Starting with this*point, "addi cronal=points om the 

inner segment of the curve were obtained keeping Z5=-1 

and sweeping along the R-axis retaining the most recent 


values =for ’pyandsA vas Siniltial approxima tons ™=ror= the 


iteratwvons. 
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CHAPTER V 


In this chapter a radial expansion which 
becomes accurate as R>0 will be shown. 

This expansion was tested in a fashion similar 
to the test procedures followed with the 
previously presented expansions. The expansion 
showed the expected behaviour in these tests. 
However it was not employed to calculate physical 
systems and there are some doubts that it can be 
efficiently used for numerical work unless an 
economical algorithm can be devised to evaluate 
two non-linear parameters in the expansion such 
that the convergence rate of the expansion is 


optimized. 
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Radial Expansion in Eigenfunctions for the 


Unasted SAtomelimit 


If R-0O the GSW-DE 
{ (\?-1)D*+2 (m+1)\D+m(m+1)-A} L(A) = 0 [401] 
for large A's becomes an Eulerian DE 


{07+ (2m+1)O0+m(m+1)-A} L(Q) ~ 0 


© =d(d/da) , |Aalz~ [402] 
whose char. polynomial has the pair of roots 
ig = +VA+% - (mt) [403] 


Based on the "dipole theorem" and the char. exponent 
€ which was determined in the previous section, the 
plus sign in [403] is the physically meaningful choice. 

We proceed to factor out a power of XA using the 
GSW-DE with coefficients as specified in [212]-[216] 

{x (x+4p)D?7+[4p (m+1)+2 (m+1-2p) x-x*]D+[ (Q/p) 


~(m+1)]x+20+ (m-2p) (m+1)-p*-A} Sw(x) = 0 [404] 
Sw(x) = x%Sn (x) [405] 


which substituted into [404] gives : 
{{x (x+4p)D?+ [4p (2a+m+l1 )+2 (a+m+1-2p) x-x*]D+[ (Q/p) 


- (atm+l) ]x+2Q0+ (m-2p) (2atmt+1)+a (atl) -p*-A}+ 
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ipalee) 
Myo eo) eee | Sn(x) = .0 [406] 


where for Sn(x) an expansion in confluent 


hypergeometric functions is desired 
Sn xe 3 gq? (n,c, 7x) [407] 


Substituting this expansion into [406] and making 


use of the relations : 


| 
5 
oO 
Ss 
+ 
Re 
Q 


xO (n,C,;x) = , yix)tlc ~2n)®(n,c,7x) 


+(n-c,)¢@(n-1,c,)7x) [408 ] 


II 


xD (n,C, 7x) n{o(n+1,c,;x)-®(n,c,7x) ] [409] 


gives the following expression 


m qn {Cn ((2(atmt1) -c, )+(Q/p)+n-a-m-1) ]x@ (nt+1,c, ;x) 
+[ (c,-2n) ((Q/p)+n-a-m-1) -n (2 (atm+1) -c, )+20+4np 
+(m-2p) (2a+m+1) +a (at+1)-p*-A]x (n,c,;x)+E (n-c,) 
((Q/p)+n-(atmt1l)) ]xd(n-1,c,;x)+04p (2atmt1—c, ) (xD) 


tdap (mita) tO (nc. >x) = 0 410] 


1 


The following substitutions are convenient at this 


point 

Gar(r) = n(atm+l+n-c,+(0/p)) (eA eled 

Gb(n) = (c,-2n) ((Q/p) tn- (atm+1) ) -n (2 (atm+1) 

-c,)+2Q+4np+ (m-2p) (2at+mt1) +a (at1)-p*-A 

[4 2a 
Gce(n) = (n-c,) ((Q/p)+n-(atm+1) ) balsa 
n= 4p (2atm+1-c, ) (a d4g 
vy = 4ap (mta) P4154 


Then making additional use of the pair of relations 
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[408] and [409], the linear terms in [410] can be 


eliminated : 


L qn { (nt1)Ga (n) ® (n+2,c);x) +0 (c,)-2(n+1) )Ga (n) +nGb (n) 
n 
tnnJ@(n+1,c);x) +L (nt1-c, )Ga(n) + (c,-2n)Gb(n) + (n-1) 
Ge (n) + (n-vn) 1@ (n,c, 7x) +f (n-c, )Gb(n) + (c,-2(n-1) )Ge (n) J 


(n-1,c,;x)+(n-1-c, )Gc(n)®(n-2,c,;x) } = fi [416] 


Collecting equal expansion terms in [416] gives the 
following five-term recurrence relation among the 


expansion coefficients in [4071] 


qn-2t (n-1)Ga (n-2) It+q [n (n-1) + (c,-2n) Ga (n-1) 


+ (n-1)Gb (n-1) ]+q,£ (n-vn) + (nt1~-c, ) Ga (n) + (c,-2n) Gb (n) 


n-1 


+ (n-1)Ge (n) J+q, ,, 0 (ntl-c,)Gb(nt1)+ (c -2n)Gc (n+1) ] 


i 


tq 4of (ntil-c))Ge (nt2) J = 0 Leola 


The choice of the parameters "a" and "cy" Sea toi erany. 
For the two limiting cases R*~ and R>0 the fastest rates 


of convergence predictably occur if 


aU ae CaS m+l1 , R>0 [418] 
a = VA+hk - (m+) , ony = 4/4A+1 +1 , R-0 [419] 


If these two parameters are chosen as in [418], 
approximating the radial w-qround=state: wi the 7 or —21 
terms, the following numbers of significant figures for 


the separation constant A are’ obtained : 
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As expected, the tables above show the same 

Chanda ceeriSeLteswas  Chesones Shown Natl O.uel ll. LO 
find the same number of significant figures for the 
separation constant as shown there, about twice as 
Many expansion terms had to be included if use was 
made of the five-term recurrence relation [417], 
although the two expansions themselves are formally 


1dentical. 


ach iuumecl TEC] mi " are chosen as in 
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[419], the number of significant figures will increase 
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This trend was empirically verified by data shown below : 


om O 20a U3 sens 04 06 0.5 an) GM A07) m0 Cee Gamal 
Q- 
5 4 4 3 3 3 2 2 2 
4 3 2 2 2 i 0 2 2 
6 5 4 4 3 3 3 3 2 
7 6 5 4 4 3 3 3 3 
; 8 6 5 5 4 4 4 3 3 


where again the radial Y-ground-state was approximated 
using) Zo-expansion™ terms, in, [4071]. 

The confluent hypergeom. expansions with three-term 
recurrence relations do not converge satisfactorily if 
RvVO. AS one can see from the table above, the drop-off 
in the number of significant figures with departure from 
the limiting case p%0, especially for small R, is so 
drastic, that it appears doubtful that the entire range, 
in which the previous confluent hypergeom. radial 
expansions did not converge satisfactorily, could be 
covered. One major problem appears to be that for small 
R we have Im(a)z0 which gives rise to an ugly singularity 


at the expansion centre 


x? = oe 2 feos (ima vin ved Stn Conan 


[420] 


A scheme with promise for improvement makes use of 


the Vaciational principle, Ussing for owathe trial stunction 
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Sw(x) % «exp (-x/2) LP (x) Vis 
where the Laguerre polynomials can be expressed 
LE (x) = @(-n,B+1;x) [422] 
For* "A", tne “approximatiom by the Rayleigh 
coefficient is equivalent to 


An i xexp (-x/2) LP (x) [x (xt4p)D?+2 (m+1) (x+2p)D 
- (3) x?-((Q/p) ~p) x+m (m+1) +20-p? Ix*exp (-x/2) LP (x) a 


i, fox? exp (-x) CLP (x) 12x [423] 


Here use has been made of the differential operator 
satisfied by Sw as shown in [404]. 


Defining the quantities , 


fox“ exp (-x/2) Lb (x) Cx (x+4p)D?+2 (m+1) 


ener a 
(x+2p) D- (4) x?+ ( (Q/p) -p) x+m (m+1)+20-p? ]x°exp (-x/2) 
LP (x) dx 
Te A ey eves Oy EO 
ee gaeal = Sy Xe EXD ( x) Ly (x) Ly (x) dx [424] 


higher orders of approximations Of | A. Canebe 
Gbtainedsin the tashion 2ndicatedvins 304) ands hs 900% 


Using the definition of the Gamma function 
fo x“exp (-x) dx = I (atl) [425] 


and writing the Laguerre polynomials as 


B BY GOT rah verdes ees) spe 
Ly (x) = on (-1) en iies A Ehal [426 ] 
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rn 
Seine AAC ((-1)7/ (4) 1) (eae) Be (ety Aca da0) 
he Te(O/ py retained) 1 ig oo Art Mere). (4 Fo 


+2 (m-2p) )+ (m-2p) (m+1)+2Q-p? IP (itj+20+1)4+(4 (j+a) (j+o+m) 


pir (itj+20) } = 0 [427] 
a = § (certyciyry (88) FZ 1)I7(4)! 
K+1,nt1  ;2, pe Ve oe Sra A 
(5) Pr (it+j+2041) [428] 


For the ground state [423] simply becomes 


pV oat [429] 


Shey silyl 
For the expansion of a given GSW-DE the rate of 

convergence is considered close to its maximum if the 

Ray lergiecOclitc1entwins U423), OrmntgneE approximations. 


are optimized with respect to the non-linear parameters 


ele aricl c,"- (*) 


Since the expansions for Sn become truncated 


for the two limiting cases R>~ and R>+0, if the 


conditions [418] and [419] are specified, the 


optimization suggested above should give for "a" and 


"Cy 


limit is approached. (x x) 


Pechescondyitions f4lsdeand) (4197s thes correspending 


(*)where the restrictions Im(a)=0 and Im(c,)=0 seem 
appropriate 
(**) Notice that for the cases [418] and L419] the radial 
expansions become truncated indicating that for the 
n-th excited state the approximation Tae eee 
is not suitable but instead higher orders of approximations 


of A (at least n-th order) seem to be required. 
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Although the algorithms mequired for thestreatment 


above are straightforward, no actual numerical 


investigation of these Speculations was undertaken. 


; aa nl 4 | A, a 7 
iwipt saene aay! - fe ose 
Jew va 


3 a 


HOE 


NEA 


1h 


5 


BIBLIOGRAPHY 


DtEratconyeu.A.7 Morse PaM., sCUlU aGaleaancerinte ir mk ane 
(1941), "Elliptic Cylinder and Spheroidal Wave 
FUNCTIONS. , J.) Wiley. Nay cand Chapman se halle London. 


SCEACCON piicA 1,0 MORSemmeD Mee meCIiT folio) gels ttl emus DeGe 
ands Corbato,. tJ. (h950) “Spheroldals Wayeurunctitonse:. 
Wiley, N.Y. 


heave iylleraas, Zelesonrietetier epiysi ka Vol. ws eels) 
Oo os 


Ree beller sete wcahlain, ee Phys cal eChenistry,.ansAdvanced 
Eueatiser, mVie a Deo ye Naw ACA em eressmclo 70) 


Jatre:G-) (2934) ) Veitsehrite fueraeny si tae ou oc 


"Handbook of Mathematical Functions with Formulas, 
Graphs, and Mathematical Tables", U.S. Dept. of 
Commerce, National Bureau of Standards, Applied 


Mathematics: series, 55, p./56, (1964). 


Erdélyi, "Higher Transcendental Functions", Mc Graw- 


HelLiaspook’ Company, snc (L953)% 


Obenheteinger,.6., bad 1, ues ebeind io DUEnOeT mod sok, 


"Table of Laplace Transforms". 


NemRry lorry Nes HOgOLInN boy oul LebineAcad mUcony 
Gening nada 29 )iep aay le to. 


MH. eMiteleman & RoE. Von Holdt, Bhys.ekhev.ep. 140, 
US ten Piao 


See eneiie ONY Ss eReEVive el 41,0000 lez te 


U2 heuer onde wl Ox, aL OOO) Phy S muetLelcmees, 
p2o47. 


ee a : a 7 

7 ae ei ‘ Ei be: dt D wt ny ve 
Vet, (ABSA CRIN Paras an ghd 7 ets | 
> eet ae] a Mae a ye 


o 


Pots eer ern. AD colt vena we +’ to: 
sacihp ate ioemiae” kee ose tae © 
Ae 


* sep hyva*) Shee" at, eae 


J my » j >t a. d 
‘ $2 on fo. eres oF ge? of a¥ ’ ‘ 
f } ’ j + ’ a) 
; y 
pP i. 4) tad Mites -® 


lier Dro tentette? hae » 


"i eo .gaead canstawt 
6). 32 4% «oP sels eee oe 


im mg 
¢ tse? vm sap" aig 


caw, coat. eames toad Lia 


2 veswd > db’ - Ube tT: omteeral oat 
aieuZineiae side tes 


ry ca 
it gina tetepiitesyeodghsee) . at fs pais: 
kT Wi afee ts /f } 


ole 


$1.9 epi) ey ~oo6gF ws 3.4)6 


Va 


APPENDIX | 


Computational Details, 
Program Listing and 
++ 
Sample Calculation for the HeH 
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Computational Details 


For a given internuclear distance the separation 
constant and the energy are varied until they match 
for the radial and the angular expansion. 


This is done as follows 


For a given R and E, (R) ands angina ala 


formulating the algebraic eigenvalue problem 
det p =A is)s=— 0 [430] 


COrresponds te looking fon thesrootssot, therciar. 
polynomial 


Wey) = [431] 
which can be done by Newton iterations 


Ag, 7 (=eats ip, (Aa) py (an [432] 


Tiel 
To expand the characteristic polynomial of a 
tridiagonal matrix is easy enough. The expansion of 
the characteristtc polynomial, Of amtridilagonal, matrix 
was previously employed to show that the Hylleraas 


and the Jaffé expansions have identical rates of 
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convergence. 
We can however equally well expand derivatives of 
the characteristic polynomial. 


For thesLirst derivative one obtains from [51a 


Bi (eaas (op, -A) Py _y (A) -B, ¥, Py_ 9 (A) -Py_y (A) [433] 
an diet tel a lelsy, Pg (A)=0 and pj (A)a-1. 


For a given set of E and R the products BLY can 


k 
be evaluated, stored and recycled in successive Newton 
iterations. 

Using Newton iterations has one outstanding 
advantage... It 1s the following : if.one, evaluates -the 
characteristic polynomial of the algebraic eigenvalue 
problem there is great danger of overflow or underflow 
due to the limited range of magnitude allowed for real 
numbers by the available software of the machine. 

However if the characteristic polynomial and its 
derivative are evaluated simultaneously one can easily 
implement a "trap" in the program which multiplies 
all recurrence terms by an integer power of 16 
depending on which one of the two problems is building 
upeselLien) since ,.nlysthe weatio Pr/Pn is of interest one 
does motuhavesto keep etrackiof allbethe multiplications 
performed sto, prevent (sdisastern 4 In practice; those 


safety features apparently require more machine operations 


than the actual Newton iterations. 
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If for the angular and the radial expansions the 
roots of the characteristic polynomial are found one 
definesea function ofp whieh vs the difference 
between the pairs or roots. Pf thisstunction 2s zero 
the desired match between the two expansions has been 
achieved. The root of the: latter scunction as) located 
using the secant method. Since the Newton iterations 
are more rapidly convergent than the secant method, 
the total number of Newton iterations is restricted 
in the’ program to not more than’ Six. These might not 
finally converge unless the secant method root search 
itself has advanced to its final converged state. 

In order to test convergence in respect to the 
number of expansion terms considered, the program 
increases temporarily the demanded numbers of 
terms of both expansions by four and proceeds to 
calculate a pair of separation constants with the 
extended expansion set using the previous final 
converged value for p. If this pair of separation 
constants differs” by more* than’ a*tolerated threshold 
SUPP ILedeiNnetheninputats Valles lseprintedsau the 
right end of the corresponding output line. 

Since the Newton iterations are very rapidly convergent, 
the algorithm described seems reasonably efficient in 


comparison with possible two-dimensional search schemes. 
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The computer program can be operated such that for 
a given energy state the separation constants and the 
energies are evaluated for many internuclear distances 
differing by a constant increment. 

If this increment is comparatively large especially 
for heteronuclear species and large internuclear 
distances, there is great danger that the search jumps 
from one state to another or flips repeatedly between 
two states. This erratic behaviour can be considerably 
Suppressed if one uses interpolation polynomials to 
predict the immediately following numbers of the 
search. 

This is implemented in the computer program using a 
weighted second order least-squares polynomial. If 
Previously at) least four points one Chesh (Re curyesnad 
converged,the new values for the separation constant 
and p are predicted giving the points increasing weight 
as one approaches the most recent point. The program 


uses the following 


Wa ey wo=1/V3 ; wi=1/V2 wel [434] 


i 4 
A weighted least-square scheme for a second 


order polynomial is 
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where x would be R and y either the energy 
(actually Es ath Ee EA) or the separation constant 
A. Since for the points considered the internuclear 
distances differ by a constant increment we can 


make the following arbitrary choice 


Se ak Sa es ek [436] 


4 
and the future point would be located at 


x = 1 [437] 


ees 2 mont 
ENE ante hha ones = ayta,ta, [438] 


is equivalent to 


eee 0.668151137905297y,-1.004453413715892y, 
~0.995546586284109y,+2.331848862094702y, 


[439] 


If less than four but at least two previous 
points were calculated and had converged, linear 
interpolation is used. 

The computer program allows the matrix 
elements or the eigenvectors to be printed ona 
separate output unit. These features are however 
NOt at abl) excitingvand| tor ithevactual usesoretche 
program at’ theybeginning Of (the sourceslisting one 
will find instructionseforelherprevaraUioneol input. 


Furthermore a sample calculation is shown 
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following the program listing. 


When using the program for numerical work one 
usually starts with a small internuclear distance for 
which the energy and the separation constant are 
approximately known by the "united-atom" limit and one 
chooses an increment for R which carries calculations into 
therueqion Ofeinterest slfirones starts: with tcor smaliean 
internuclear distance one occasionally encounters 
convergence problems for the radial expansion. 

The computer program is rather fast. When for a great 
number of internuclear distances energy calculations were 
performed for the Hen (6,271 state, approximately 58 
milliseconds were needed per point on the E,(R) curve. 
The machine employed was located at the University of 
Alberta which presently (1977) owns a hybrid consisting 


of IBM 370 software and an Amdahl 470 CPU. 


Referring to the same calculation with the same 
machine the following histogram of distribution of CPU 
time among the main program and its subroutines was 


recorded : 


MAIN (Main program) OSs 


FDIFF (Function defined as the difference 
between the separation constants 
for the radial and the angular 


expansion for a given energy) 0.978% 
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DETA 


DETHET 


DETR 


DETR 


DETER 


EIR alt 


INTPOL 


INTPL 


DSORT 


(Galculation ofthe matrix 
elements for the homonuclear 


angular expansion) 


(Calculation of the matrix 
elements for the heteronuclear 


angular expansion) 


(Calculation oF the matrix 
elements for the radial 
expansion) 

(Newton iterations) 
(Output of tridiagonal 


matrices) 


(Calculation and output 


of eigenvectors) 


(Controle subroutine. fot 


interpolation) 
(Interpolation subroutine) 


(Double precision IBM systems 


Subroutine for the square root) 
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0.000% 


12 ae 


e035 


(ei eVAW: 


0.000% 


0.000% 


0.030% 


020712 


0.438% 


The corresponding object program was prepared 


by the “FORTH” compiler specifying unmestricted 


OpbiIMi Zadtwon OptLon. 


The expansions used in the computer program are 


Homonuclear angular expansion 


(Generalized Legendre polynomials) 


Formula [14] 
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Heteronuclear angular expansion 


(Ultraspherical polynomials) 
Formula [38] 


Radial expansion 


(Hylleraas expansion) 
FOrmuULamy (3) 


The formulas indicated refer to the recurrence 
relations for the corresponding expansions. 

The matrix elements for the homonuclear angular 
expansion are as indicated in [14] i.e. A, is 


dividedibyepc. 
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Sig VA Bee THIS PROGRAM PERFORMS CALCULATIONS FOR 
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METHOD ; 


PEPER SENN, MASTERS LHESLS CU977) , UNOVERSETY OF ALBERTA, 


CANADA. 


LEME. S 
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PEE Wee see a OUlP UR ON CURE CNter ome s CCROU LED SEN 
BLOCKS WHICH ARE SEPARATED BY HORIZONTAL SPACINGS. 
VES SDE RM UNE S! Lop Sly OR CHEST BLOCKS ANT Is 
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PROGRAM @ Diy NORM AU isc Dao MOR MAGE GS) is 


LOCOMDT coe e SOR He HOMONUCTL EAR CAS i ten UT Se tea) 
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FOR “LOGODD", ANTIBONDING STATES WILL BE CALCULATED. 


Neves es hs DLMEN STON OPV RRE MARA eB 
CONSIDERED FOR -PHE ANGULAR EXPANSION. EF “'N! IS 
ZHRO OR NEGATIVE THE CORRESPONDING CARD IES THE 
LAST ONE. THE; PROGRAM SETS N = { Nt. 


DIR Mieve are Hah, SAME AS CON MOSS BU Te SRB RUB ROR teres OS Eine: 
RADIAL EXPANSION. 


M .«<ee. THE ANGULAR QUANTUM NUMBER AROUND THE 
INTERNUCLEAR AXIS. 


NRAD «secs. LHE NUMBER OF INTERNUCL RAR DISTANCES 

WITH AN INCREMENT “RDIFF" FOR WHICH ENERGY CALCULATIONS 
ARE DESIRED. IF "NRAD" TS ZERO OR NEGATIVE THER PAIR 

OF MATRICES WITH THE FINAL CONVERGED SEPARATION 
CONSTANTS WILL BR PRINTED ON OUTPUT UNI "7" IN 
HEXADECIMAL CODE. 

THE DIMENSION OF THE OUTPUT MATRICES WILL BE N. 

THE PROGRAM SETS NRAD={NRAD]. 


MAXIT ..... THE MAXIMUM NUMBER OF FUNCTION EVALUATIONS — 


PERMITTED TO LOCATE THE ROOT WITHIN THE ERROR BOUND 
SEC Gren by eG pat. 

EE MAMTIT ©S CERO OR NECATIVE THE ERGENVECTORS OR 

THE FINAL CONVERGED ENERGY AND SEPARATION CONSTANT 
VOLE BE CALCULADEO AND PRINTED ON OULEUT UNL TSie 
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By selec) DHE GNIDTATL CHONCE FOR THE RERCTRONICG BNERGY. A 58 
"E@ HAS TO BE LESS THAN ZERO. es) 
LY 360 

R eeeee INITIAL INTERNUCLEAR DISTANCE FOR WHICH THE A 61 
ENERGY AND THE SEPARATION CONSTANT WILL BE CALCULATED. A 62 
(SEE NRAD £11!) yA UeyS} 
A 64 

ROLE cc. LUN THERNGCL EAR DUSTANCE TUNG REM ie (CS A 65 
NRAD !!4) a 6S 
A 67 

A aceee INITIAL CHOICE FOR THE SEPARATION CONSTANT A 68 
POR WHECH THA CALCULATIONS ARE EXPECTED 80 CONVERGE iN sey 
TO THE DESERED ENERGY STATE. A 70 
. el 

COMMS sic aie UNC RE WENT PORE tet CO SAN) ere ROG A te 
FINDING PROCEDURE USING A MODIFIED SECANT METHOD, A PS 
WHERE C = -0.5*2*R*R. ‘ A 74 
f A 165) 

UMNONM KA Go THOR TIAP INDY GnRONONS Ole UERUE IRIE, (Oy meen A 716 
ENERGY. h ae 
A 72 

a a a = = = = = - PAS, 
4A ec 

THE NEXT CARD CONTAINS A LOGECAL CONSTANT AND THE A 81 
NUCLEAR CHARGES. THE LOGICAL CONSTANT IS SE a A 82 
HOMONUCELEAR SPECIES 2S BERING CALCULATED AD) so aes. A 83 
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A 85 

a ee eee A a6 
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A B88 

THE LNPULT TS GREANTSED AS FOLLOWS <= A 84 
A YG 
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K 932 

THES DATA IS READ FROM INPUT UNIT "5S" UNTIL “NY TS READ A 93 
EN AS ZERO OR NEGATIVE. A ou 
95 

Tiga JMO AUS NDAD) BIBS A G6 
A 97 

FORMAT (14, 14,13,127215,5810.0,2810. 3/01, PS. 0, 10. 0) A 98 
A 99 

A 100 

Soe ce PETER SENN , UNIVERSE TY OF ALBERWA SS SoS AOLA 

A 102 

a AUGUST 197 SS eae 

A 70u 

A 105 

IMPLICIT REAL*8 ee =i) pee N eh GE Ree 2 le IN}) A 106 
REAL*8 UP1(200) ,DIAG1(201), LOW 11200), Ah 107 
# UPZCZ0Ci D1AG2 (201) » bow? ( 00) A 108 
LOGICAL*1 LOCODD, LOGS, LOGU, LOGIT, LCGVEC, LTR,LSTOP, IFLAG, Tee) A 109 
COMMON. /DIFF/ DIAG1, DIAG2, UP, UPS, LOM LONG | Rh pall 5 Acz pide al only polices pe Sgt: 
#/EVODD/LOGODD/COUN it /ISL/STOP/IFLAG CHARGE/ZA, ZB Mh 144 
ee (A1,A) eae 


GATR 10,11,13, 18, 16, 17/0, lip3e4, Oy 7/7; OND, HALE/ 1.0005 0. SD0/, LSTOR/. watts 
APD esta R/. TRUE. /, CCUT/Z3A10000000000000/ A 174 
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READY (5S teh TLING A a%5 
Cote=e TO reer. 
WRITE (6,14) Ao dag 
READ (5,13) LOGODD,N,NR,M,NRAD, MAXIT,E,R,RDIFF,A,CDIFS, CUTOFF, A 118 
#L22,ZA,2B Nh We 
CUTUFF = 0.125D0*CUTOFF A 120 
fs) OR A 121 
Sips ye O00 a oo 
CT = ONE+CUTOFF A 123 
TD (HG IIR TE) SOR = MRT A 124 
IF{LSTOP)N = -N Ades 
LOGVEC.= sFALSE. A 126 
LP (MAXIT. BL. 21 LOCVEC- = a TRUE. A 127 
IF(LOGVEC) MAXIT = -MAXIT A126 
LOGU = .TRUE. f KR 129 
LF (NRAD.LT.L1LOGU = . FALSE. A 130 
IF(.NOT.LOGU) NRAD = -NRAD A t24 
Ni = N-121 A 432 
N4 = N+T4 j A 133- 
N4M1 = NtI3 A 134 
NIR = NR-I1 A185 
N4R = NR+T4 A 136 
N4M1R = NRtI3 A 137 
Resa KI A 138 
C = -HALF*E*R*R A 139 
DOr ti I=h 1p N RAD A 140 
SJ = SJ+ONE A 141 
LOGS = sFALSE. A 142 
LOGIT) =i1,0GS A 143 
TY Se A Yua 
FUNG = PDIFE (XN eth) A 145 
CDIFF = 8.0DO*DABS(FXN) A 146 
CDIFF = DMINA(CDIFF, DABS(CDIFS)) A 147 
TPCODSPr. Leer) enarr = CpaLEs A 148 
XNM1 = XN+CDIFF A 149 
FXNM1 = FDIFF(XNM1,17) A 150 
FXNP1 = FDIFF(XN-CDIFF,I7) Sa 
IF (DABS ( FXNP1).GT.DABS(FXNM1)) GO TO 2 AL Sez 
FXNM1 = FXNP1 A 153 
XNM1 = XN-CDIFF Ay Se 
Lean a0 K 155 
[FLAG ="LOGS A 156 
BOW) JVR=11;, MAS IT be 
[Ff (IFES) Go Te 1 A. 158 
IF ((XN.EQ.XNM1).OR. (FXN.EQ.FXNM1)) GO TO 4 A 159 
XNP1 = XN-PXN*((XN-XNM1)/( FXN-FXNMN1)) A 160 
XNP1 = DMAX1(XNP1,CUTOFF) A 161 
Go TOUS A 162 
XNM1 = XNM1¥* (CT) A 1603 
FXNM1 = FDIFF(XNM1,17) A 164 
ISK = ISK+tI1 A 165 
IF (ISK.LE.MAXIT) GO TO 3 A466 
LOGS =). LRUE. i SY 7 
corre 8 A 168 
FXNP1 = FDIF?P(XNP1,15) A 169 
XDIFFS=XNP1-XN AO170 
XNM1 = XN BAG 
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FXNM1 = FXN 

XN = XNP1 

FXN = FXNP1 

IF (LOGIT) GO TO 6 

IF (DABS (FXNP 1/A) #DABS ((XN-XNE1)/XN) WET. CUDUFF) LOGIT<=) .TRUES 
GO TO 7 


IF (DABS (FXNP1/A) +DABS ( (XN-XNM1) /XN) .G2.CUTOFF) LOGIT = . FALSE. 


Die (OCI AiGOm Lee 

IE CURE ARMOR TU RIS UEONE;S) = ara 

Cea kN 

IF (DABS (FXNM1) ~.LT.DABS (FXN))C = XNM1 

CHECK ON FINAL CONVERGENCE. 

EP {L2Z2) CALL) DEAS (UPA DIAG, LONA, 6, Na) 

WE {G8 WON AVA mATy  CUMsi, IDS URE ND (UP1,DIAG1,LOW1,R,C,M,N4) 
CALL DETR (UP2,DEAG2, LOW2 ,R,C, M, NER) 

BORG IR 17 ae 


a1) eee) AY = —-DETRE (UP1, DIAGT, LOW1,-Al1/C,N4, Nan 1, L&R) *Cc 
TPF (.NOTlLLZAZ) Ad = DETRI(UPT,DEACI, LOWT ANG, NOM, LER) 
A2 = DETRI(UP2, DIAG2,LOW2,A2,N4R, NMR 7p le) 


IOS ES GROIN S 

TD he ae RC Us 

WN = DABS((A1-A2) /A1) 

LCN. Cr. COLOR) LOGS = PRUE. 

Ish = oe 

EF (1.06 <8. 2081) GO TO 10 

OUTPUT OF HE sil Aa Gms). 

IF (LZZ) CALL DE RA (We 1, DIAG LOM Gy i, NY 

Lt toNOl. L247) CALL DET HED UR 1 PDA G Tea OW dan, © lp) 
CALE DRLR CUR2, DEAGZ 


GA, LOND, hi, C, a, Nit) 
PEGA CALL SUBD ER (CO Odea AG Ae TON cane Git pele 1) 


c 
Le Ca NOS sc) CAL iD Drees (ULE. 1, Daten eens A1,N,N1) 
CALL DETER (UP2,DLAGZ,LOWZ, A2, NK, N1R) 
AT = HALF*(A%+A2) 
SG be) en) 


ETOT = EtONE/2 

TE ¢.NOT.LOGS) WRITE (6, 15) ISL, 2; My, 
IF (LOGS) WRITE (6,15) ISL,R,M,ZA,ZB,5 
DATA ICOUNT/O/ 

ICOUNT = ICOUNT+I1 


| WP ae Als) Palin is EL O LIE ad VAY 
7 SOT FA WN 


IF (C((2COUNT/ILINE) *LLINE).EQ.TCOUNT) -AND. ((J.NE.NRAD) .OR. (.NOT. 


#LSTOP})) WRITE (6,16) 
IF (J.LT.NRAD) CALL INTPOL (A1,C,R¥*R, (RtRDIFF) **(2),5,LOGS) 
hee = At 
Lh ((.NOT.;LOGVEC) .OR. LOGS) CO TC 11 
CALCULATION AND OUTPUT OF THE EXPANSION COBFFICIENTS. 
LF ELS}. CALE DETR (UP 1, DIAGT, LOWT,C, 4") 
IF (L2Z) CALL TRI (LOW1,DIAG1,UP1,DIAG2,-A1/C,N) 
IF (-NOT.LZZ) CALL DETHET (UP1,DIAG1,LOW1,R,C,M,} 
IF (.NOT.LZZ) CALL TRI (LOW1,DIAG1,UP1,DIAGZ,A1,N 
CALL DETR (UP2,DIAG2,LOW2,R,C, M, NR) 
CALL TRI (LOW2,DIAG2,UP2,DIAG1,A2,NR) 
R = RS+SJ*RDIFF 
IF (.NOT.LSTOP) GO TO 1 
WRITE (6,17) 
STOP 
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FORMAT 
FORMAT 
FORMAT 


C25} 
(fp Oe, 3,02, 2b oS P10. 0,8 10.3845 1 Pos 0, P10. 05 
(TTA =, 20x paae eee CENTER, (X%, I2ZHONE=EL POURON, 1x, 


#7HPROBLEM/1H-,8X,1H+,5(19(1H-) , TH+) ,/9X, 1H ,5 (19%, 1H1) ,/9X~ WH] 2k, 
41 OHINTERNUCL. »1 4, SHDIST., 1X, 1849 3X, 3HPHI, 1X, 7HOUANTUM, 1X, 1H#,3X, 


#1H] , 3x 
#1X, 108 


, OHELECT., 1X, OHENERGY, 3X,1H{, 4X, 5HTOTAL, 1X, 6HENERGY, 3X, 181, 
SEPARATION, 1X, GHCONST., 1X, THIS /O RINE eo C19 Le 10H 


#/9X, HL, 19X, 1H |, 3%, SHNUCL. , 1K, THCHARGES, 3X, 1H{,-3(19X, 181) -/9X- 
#1H1,5(19X, THI 72x, 1H+,5 (19 (1H-} H+) ,/9%, 181, 5 (19%, 1H) ) 


FORMAT 


FORMAT 
#)) 

FORMAT 

END 


THIS FUNCTION EVALUATES TH 
TWO EIGENVALUES OF INTERES 


os 3X, VEL, P 17. sie geek 1H pI4,2P6.2,3%, 1H ,3(P17.12, 2X, 


i UE C2 ee Oe 3 OX SBAE PE) 


(9X, 1Hf,5( 19%, 1Hf} ,/9X, 1H+,5 (19 (tH-) , TH+) ,/9X, 1H], 5 ¢19X,18] 


(9X, 1H{,5(19X, 1Hf) ,/9X, 1H+,5(19 (1H-) , 1H+) ,/1H1) 


TFFERENCE BETWEEN. THE 
THE PAIR OF 


1S 
ey 


TRIDLAGONAL MATRICES 


REAL FUNCTION FDIFF*8 (C,NEND) 


TL eB IGA EAE 


REAL*8 


REAL*8 (A-H,O-Z) , INTEGER*2 (I-N) 
UP1 (200) , DTAG1 (207) , LOW 1 (200), 
UP2 (200) , DIAG2 (201) , LOW2(200} 


“OMMON /DIFF/ DIAG1, DIAG2,UP1, UP2,LOW1,LOW2,R,A1,A2,M,N, NR, LZZ 
LOGICAL*1 Sy Cee: 

EQUIVALENCE (A1PREV,A2PREV) 

DATA CPREV,RPREV/2 ee 123459876D71/,11/1/, PASS/1.0D-14/ 


LOGA = 
LOGO = 
pie ue 
LOGO = 
CPREV 
RPREV 


sig gEwE 


Tet. NOL. 


«FALSE. 
LOGA 
HBO. PREV AND (REO RP REV) i CO TOs | 
> TRUE. : 
eS 
R 


Zy) CALL DETA (UP1,DIAG1,LOW1,C,M,N) 


L224) CAL. DERE T CUP 1D Ac 1,108 1p R-C, M,N) 


CALL DETR (UP2,DIAG2,LOW2,R,C, M, NR) 


N1 = 
i= 
TP (he 


als 
NR-11 
Z) Al —-NETRI(UP1, DLIAG1,LOW1,-A1/C,N,N1, LOGO) *C 


IF (.NOT.LZZ) Al = DETRI(UP1,DZAG1,LOW1,A1,N,N1, LOGO) 
A2 = DETRI(UP2,DIAG2,LOW2,A2,NR,N1K,LOGO) 


ATPREV 
LOGO = 
LOGA = 
DOZ. L 


= Al 
ah Ada oiine 
LOGO 


=T1,NEND 
EF (22) Al = =DETR 
Tee NOW. Zia) Ante 

IF (LOGA.AND. (DABS¢( 
IF (DABS ((A1-A1PREYV) /A 


te ere ea 
TRI (UP 1, DIAG1, LOW1,A1,N,N1, LOGO) 

oer tite at GO To 3 
1) 


~uc. PASS) LOCA = se TRUE. 


AVPREV = At 

A2ZPREV = A2 

LOGA = Sigs dec 

Ome my =e tite helt 

ea DRG (UP TAG2,LOW2,A2,NR, N1R,LOGO) 


TF (LOGA.AND. AROS eS A2PREV) /A2) «LT. PASS) ) GO TO 5 
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IF (DABS ((A2-A2PREV) /A2). 


A2PREV = A2 


Lf.PASS)LCGA = .~TRU 


TDA SSS Au ae 

RETURN 

END 

GETS) SUBROUTINE GENERATES THE MATRIX 

ELEMENTS OF THE HOMO NUCL LEA R ANGULAR EXPANSION. 
SUBROUTINE DETA (UPPER,DIAGON, LOWER,C,M,N) 


IMPLICIT REAL*8 
REAL*8 DIAGON (2) 
INTEGER*2 


M,N 


(AH, O=2)s 7 
, UPPER (2), LOWER (2) 


INTEGER*4 (1-N) 


LOGICAL*1 LCGICoO, LOCIEI 


COMMON /EVODD/ LOGICI 
DATA ahi Ata8/ 

FUP({LPM1I1) = (DFLOAT 
FLOW(LMMI1) = (DFLOAT 
FDEAG(LI1) = ({(DFLOAT 


# (LM M*LPM) /DFLOAT ¢LL1¥LLM1) 


LOGuUGO) = 
L= MM F 
TEPC 
#LOGICT) ) 
Ie S04 
COMO 

OG dC Om eh Al Sie 
UPPER LOC) = 

DIAGON (LOC) = 


- PRUE. 


(M/2}) . BO. fl) 
Te 


£ = L+i2 

EOC = LOCFIE1 
LPM = Lt+M 

LO = Lah 

ily we weloreattate| 
is = ilk 
LPM = iP Mis 
ESN ee ie 


Ee (LOGICOj GO TA, 1 
DEAGON(LOC) = 
LOWER(LOC-1I1) = 
IF (LOC.EQ.N) RETURN 
UPPER (LOC) 
GO T0,.2 
END 


(LPM 1# (LE 
(LM Ms 


(iM N11) ) /DELOAT (L 
((MM+21) 
Jas 


LPM1) /DFL 
DPLOAT (i 


«AND, (. NOT. LOGICT) ) . OR 


FUP (LPM1+21) 
FDIAG (L+11) 


FDIAG(L+I1) 
PLOW (L 


MM-T1) 


= PUP (LPM1+I1) 


PM11I1)) /DFLOAT (LL3* ( 
LM Deeg 
OAT (LL3*L 
EEA) ZAC 


#(((2* 


(LL3t+E. 


(N72)) - 
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THIS SUBROUTINE GENERATES THE MATRIX ELEMENTS D 3) 
OF THE HETERONUCLEAR ANGULAR EXPANSIONS. D 4 

D 5 

SUBROUTINE DETHET (UPPER, DIAGON, LOWER, RST, P2,M,NDIM) D 6 
IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*4 (I-N) D 7 
REAL*8 UPPER (2) ,DIAGON (2) , LOWER (2) D 8 
INTEGER*2 M,NDIM D ) 
LOGICAL)  MOGTEO, TRLAG D 10 
COMMON /STOP/ILFLAG/CHARGE/ZA,2ZB Dent 
RADA F414 D112 
FUP(F1) = F1*(B+PT2*DFLOAT (N+M22)) Des 
FLOW (F2) = F2*(B-DFLOAT(N11) *PT2) Dee 
LOCTCOR=se TRUE D 4S 
PRS (CP! L P.O. 0004 Go T1613 . Dieedio 

R = RST*O.5D0* (ZA-ZB) Dae 
MNT = M* (M421) : D 18 
M2 = M+M : D 19: 
M21 = M2+71 db) Ae 
M22 = M274251 D 25) 
M23 = M22471 DR 22 
M241 = M2-11 D238 
P= HSORM(P2) i) QE 
PT2 = pep fp 25 

B = R+R-P*DFLOAT (M22) D 26 
eee ig i Dy 27 
rit = D 28 
GOMEO) 2 : D 29 
TOGCLCO.= SPRL Sm. ; D1 30 
UPPER{(N1) = FUP(DFLOAT (N+M21) /DFLOAT(N2+M23) } Dee Si 
DIAGON(N1) = DFLOAT (MM14N* (N4+N%21))-P2 Dp 32 
Ze =e IN De ss 
N = N1 D 34 
ia) =e eine ) AS 
N2 = N+tN D 36 
iF <LOGICO) GOTO. 1 De eid 
DIAGON(N1) = DFLOAT(MM14N* (N4#M21))-P2 D 38 
LOWER(N) = FLOW (DFLOAT (N) /DFLOAT (N2+M2M1) ) D 39 
IF (N.EQ.NDIM) GO TO 4 Dd 40 
UPPER (N1) = FUP (DFLOAT (N+M21) /DFLOAT (N2+N23) ) D 47 
Ele) 40) D 42 
Daan (oreo) D 43 
ITFLAG = LOGICO D 44 

4 RETURN D 45 
D U6 

5 FORMAT (1H-,9X,6HSEARCH, 1X, 5SHABREA, 1X, 3HFOR, SHPOSITIVE, 1X, D7 
#8HENERGIES, 1X,9HSPECIFIED, 1X,3H!!!,/1H0) D 48 
END D 49 
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THIS SUBROUTINE GENERATES THE MATRIX 
ELEMENTS OF THE RADIAL EXPANSION. 


SUBROUTINE DETR (UPPER, DIACON, LOWER,RST,C,M,NDIM) 
IMPLICIT REAL*8 (A-H,O-Z) , INTEGER *Y (I-N) 
REAL*8 DIAGON(2), UPPER (2), LOWER (2) 

INTEGER*2 M,NDIM 

LOGICAL#] 10GICO, 1 FLAG 

COMMON /STOP/IFLAG/CHARGE/ZA,2B 

DATA I1/1/ 


FUP(NNM) = DFLOAT (NNM) -DFLOAT(N) *SQRTR 

FDIAG(NMM1) = DFLOAT(NM4+N1) *(SQRTRC- DEL aE aR Ce RET EMT 
SQC¥*DFLOAT (2* (NM+N1))) 

FLOW(NIN1IH) = DFLOAT(N1IN1M) -DFLOAT(N1-M) *SORTRC 


Lecrce = , TRUE. 

GE (CL LT.0..0D0} €o To 4 
R = RST*0.5D0* (ZA+ZB) 
SOC = DSORT(C) 

SQRTRC = R/SOC 

RRC = R+R-C 


M1 = M+T4 
N= M1 
Toei =a 1 
GO TO 3 


LOGECG = SEAUSE. 

UPPER (LOC) = FUP(N*NM) 
DIAGON (LOC) = FDIAG(N1*M1) 
N = Nil 

LOG = LOcti} 

Slee it 

NM = N-? 

eS CLOGHe oO) mGO™ LO 1 
DIAGON(LOC} = FDIAG(N1*M1) 
LOWER (LOC-11) = PLOW (1 (N1—M) } 
iP (LOC.BO. Nps Mh) GO. TO 5 
UPPER (LOC) = FUP (NK *NM) 

GO Tor 2 

WRITE (6,6) 

IFLAG = LOGICO 

RETURN 


BORMAT “(lH-,9% OHS BARCH, UX, SHR EREA x, SHnOR,, 1X, on PO oueii by, UX, 
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THIS FUNCTION PERFORMS ONE NEWTON ITERATION TO 
FIND ONE ROOT OF THE EQUATION DET{A)=0 , WHERE 
"A" TS A TRIDIAGONAL MATRIX. 


LOGIC} 


REAL FUNCTION DETRI*8 (B,A,G,Z,N,Ni,1 
R*2 (I-N) 


PMP CRCET REALS. (k=, 0-7) 9, UNTEGH 
REAL*8 A(2),B(2),G(2),LOWLIM 
LOGICAL*1 LOGIC 

DATA IONE/1/,CUT/2Z3310000000000000/ ,UPLIM/25519000000000000/, 
#LOWLIM/Z2D10000000000000/, FACTOR/Z3 71 0000000000000/7 

RES (NOT. LAG UC) iGO To 2 

DO 1 J=1ONE,N1 

BAI Wa= 8 (dt ac. () 

DRS = TO @ 

DETRI = PR2 

DERIVA = PR2 

PR1 = A (IONE) -2 

PRD Oat ee Pe 


PRDOT2 = 0.0D0 
IF (N.FQ.IONE}) GO TO 9 
J = IONE 


J = JtIONE 

AO = A(J) -2 

Al = B(J-ION®) 

OVERFLOW PROTECTION SECTION : j 
IF ((DABS(DETRI) +DABS(DERIVA)).LT.UPLIM) GO TO 4 

IF ( (DABS (2R1) .LT. LOWLIM) «OR. (DABS (PR2). LT. LOWLIM) OR. 


# (DABS (PRDOT1).LT. LOWLIM)-OR. (DABS(PRDOT2) ~LT.LOWLIN)) GO TO 4 


PR1 = PRI*FACTOR 

PR2 = PR2*FACTOR 

PRDOT1 = PRDOTI* FACTOR 

PRDOZZ = PRDOT2Z*EACTOR 

EO TOs 7 

IF ((DABS(DETRI) «GT. LOWLIM) .AND. (DABS (DERIVA) .GT.LOWLIM)} GO 
IF Cea gy ek Rd ere LL slg ces il oan bee 1 OR. 
# (DABS (PRDOT1) .GT.UPLIM).OR. (DABS (PRDOT2) .GT.UPLIM)) GO TO 7 
PR1 = PRI/FACTOR 

PAs PR2Z/PACDOR 

PRDOT1 = PRDOTI/FACTOR 

PRDOT2 = PRDOT2/FACTOR 

DPE A=" Ose Ris ee Re 

DERIVA = AO*PRDOT1- epeeeeee 2=PR | 

Ee CUE OoN sO. (DABS (2+ DE RAVAN ADE dines) aint Um} CO Ole: 
PR2 = PR1 

PRDOT2 = PRDOT1 

PR1 = DETRE 

PRDOT1 = DERIVA 


GO TO 3 
DETRI = Z-DETRI/DERIVA 
RETURN 


DETREI = PRI 
DERIVA = PRDOT1 
GO TO 8 

END 
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THIS SUBROUTINE PRINTS A TRIDEAGONAL MATRIX 
IN ITS HEXADECIMAL FORM ON I/O UNIT "7", 


SUBROULINE DETER (A,D,B, £i6, N, N11) 
IMPLICIT REAL*8 (A-H,O-Z) , INTEGER¥2 (I-N) 
REALS “A(2) 7 DiC2 ry BE2) 

COMMON TAGODNG 7 sakS i: 

DATA Lt717 

WRITE Ef ISL 

WRITE: (7,2) ETC 

DOM Smt ani 

WRITE? (7 72y Sha) 0.0) 78 (4) 

WRITE (7,3) DW) 

RETURN 


FORMAT (1X%,3216) 

ORMAT (17X%,216) 
FORMAT (1X,15) 

END 
THIS SUBROUTINE EVALUATES THE EIGENVECTOR OF A 
TRIDIAGONAL MATRIX BELON 2k TO A. SPECIFIED 
EIGENVALUE. THEN THE EIGENVECTOR IS WRITTEN , 
ON THE I/O UNIT "8" -IN Save IMAL CODE. 


LOGICAL*1 LULY 

COMMON /COUNT/ ISL 

DATA <01,92/1,27, ONB/1.0D0/ 
TsO) p= 8 te ASE 

DO 1 J=Ii,N 

DING (J) = DIAG (J) -PIG 

B(I1) = ONE 

SU TSAO AUN GO INO), 

VOW K=T 2 Nn 

IK = K-if1 

RATIO = -SUB (LK) /DIAG (IK) 
DIAG (K) = DIAG{K) +RATIO*SUP (IK) 
DOW 3 k= i257. 

i aR IG, 

B(IK) = - (DIAG (K) *B(K))/SUP{(K) 
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INS 


DO 6 K=12,N 

NUMBER = K 

if (NUMBER. EO. N) GO TO 6 

IF (DABS (B{K+I1)})-DABS(B(K))) 4,4,5 
UU RU. 

COME GMG 

EPP CLULU GeO; 20 (7 

B(I1) = B(I1)+B(K) *B(K) 

B(11) = ONE/DSORT(B(I1)) 

DO 8 K=12, NUMBER 

Bh) =) B CK) Bi Et) 
LPN. HOL EL): NUMBER = 
WRITE (8,12) ISL,NUMBER 
DO 10° K=I1, NUMBER 

WRAL Sera lo) BAR) 
RETURN 


1 


FORMAT (1X,2156) 
PORMADY (xX, 225) 
END 


SUBROUTINE TO CONTROL INTERPOLATION SECTION. 


SUBROUTINE INTPOL (A,C,RR,RD,J,LOGS) 
IMPLICIT REAL*8 (A-H,O~-Z) , INTEGER*2 (I-N) 
REAL*8 APREV (4) ,CPREV (4) 

LOGICAL*1 LOGS,LG, LAY 


DATARS Op Eg et JOR lent / wD RN iC PR ee Ore ONO 


iP Suocs), co ke 2 


LAY = .FALSEY 

JS = JS+T1 

EPG Ook 1) ase = 1 
PFA(S. Geli LAY = seRure 


ie, = . hese. 
LEQSS sit. 14) Lee = en ine 
APREV(I4) = A 

CPREV({I4) = C/RR 

CALL INTPL (APREV,LG,LAY) 

CALL INTPL (CPREV,LG,LAY} 

IF (LG.AND. (.NOT.EAY)) GO TO 3 
A = APREV(T4) 


GC = (GCPREV (24) RD 
RETURN 

JS = 10 

Gea CRD RR 
GOO on 

END 
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SUBROUTINE TO FIND THE INTERPOLATED INITIAL VALUES. 


SUBROUTINE INTPL (Y¥,LG,LAY) 

IMPLICIT REAL*8 (A-H,O-Z) , INTEGER*2 (I-N) 

REAL¥*8 A(3),Y(4) 

LOGTCAL® 1 LG Lay 

DATA A/-1.004453413715894DC,-9 .955465862841105D-1, 


AG Sis} SON SKsie 2 (CSIN F/O) W077 Ib 5 A ah sh abl 24) eh YZ 


iy (GWE) EC) ne z 

NOME SS (gots WS) 1 Shi GU ay2S) 7/S) oy {lye (LAE 4) 
DOMINI= 2704 ; 

YNEW = YNEW+A(J-71) *Y (J) 

DORS BIT ens. 

Y(J) = Y(J+t+T1) 

eI? (ie) (EG) ao) S 

Y (04) = YNEW 

RETURN 

SEIS (UIA NOTRE) SS NCEE SH) se VETCIE SH) Nea) 
GO TO 4u 

END 
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following organization : 
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APPENDIX I] 


Separation Constants for the 


Angular and Radial Y-states 
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APPENDIX IT] 


Determinant of a Pentadiagonal Matrix 


For a pentadiagonal NxN matrix the following 


substitutions are made 
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Then elimination of the elements Miao ; can be 
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achieved by the following algorithm : 


(1) penn): 
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Then we have what could be called an upper Hessenberg 
band matrix whose determinant is readily expanded. 
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where D is the determinant of M. 
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However the matrix could be brought into triangular 


form. The determinant of such a matrix is 
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This is achieved by the following algorithm 


n=0 , D=a 
(1) n=n+1 
if n=(N-l) go to” (2) 
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This algorithm requires essentially 8N multiplications 
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whereas the previous algorithm required 10N 
multiplications to evaluate the determinant. 
Both algorithms were coded in Fortran IV as 
shown below but no extensive comparison of their 
performance was made. 
The arrays in the subroutines contain the 


101 lowings: 


E(1) some eb (bee AL (i) =a 89 
BE(i)=8, , GA(i)=y, iS Gal 


and NM1=N-1, NM2=N-2. 


The subroutines shown are for complex algorithms 


and actually evaluate the norm of the determinant. 
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REAL FUNCTION PENT*8(C,B,AL,BE,GA, NM2) 
IMPLICIT COMPLEX*16 (A-H,O-Z) 
COMPLEX*16 C(2),B(2) ,-AL(2), BE(2),GA(2) 
REAL*8 CDABS 

DATA I1/1/ 

J=0 

DM2=(1.0D0,0.0D0) 

DM1=AL (11) 

D=AL (11) ¥AL (2) -B(1I1) *BE(I1) 

J=J+T1 

J1=J+1I1 

J2=J1+T1 

SIGMA=C(J) /B{J) 

B(J1)=B(J1) ~-AL(J1) *SIGMA 

AL (J2) =AL (J 2) -BE(J 1) *SIGMA 

P5=D*AL(J2) -DM1*B(J1) *BE(J1) +DM2*GA (J) *B(J) *B(J1) 
TF (J. 2O.NMo) (GORTOR2 

BE(J 2) =BE(J2) >GA(J1) *SIGMA 

DM2=DM1 

DM1=D 

D=P5 

GO TO 1 

PENT=CDABS (P5) 

RETU RN 

END 


REAL FUNCTION PENT*8(C,B,AL,BE,GA,NM1) 
IMPLICIT COMPLEX*16 (A-H,O-2Z) 
COMPLEX*16 C(2),5(2) ,AL (2), BE(2) ,GA (2) 
REAL*8 CDABS 

LOGICAL*1 L1 

DATA I1/1/ 

L1=. TRUE. 

N=0 

PS=AL({1I1) 

N=N+11 

IF «(N.EQ.NM1) L1=.FALSE. 
N1=N+tI1 

PE eCeNOU.L1)eCOMlOn2 
N2=N1+TI1 

SIGMA=C(N) /BN) 

B (N1)=B(N1) ~SIGMA*AL(N1) 
AL(N2) =AL (N2) ~SIGMA*BE(N1) 
BE(N2) =BE(N2) -SIGMA*GA (N1) 
SIGMA=B(N) /AL (Q) 

AL(N1) =AL(N1) -SIGMA*BE(N) 
BE(N1) =BE(N1) -SIGMA*GA(N) 
P5=P5*AL (N1) 

TRACE) GOmTOn 
PENT=CDABS (P5) 

RETURN 

END 
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"The only true wisdom", the caribou 
wizard told the explorer Rasmussen, 
"lives far from mankind, out in the 
great loneliness, and it can be 
reached only through suffering. 
Privation and “suffering alone can 
open the mind of man to all that is 
hidden to others". 


H. Osterman, Report of the 
Firth, Thule 
Expedition 
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